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TECHNICAL  REPORT  SUMMARY 


The  objective  of  this  project  is  to  combine  a  number  of  recent 
advances  in  finite  element  theory  and  computer  technology  for  analyzing 
cavities  and  structures  in  rock.  This  computer  program  applies  to  general 
three-dimensional  structures,  considers  nonlinear  material  properties, 
time  dependent  properties,  gravity  loading  and  sequence  of  excavation  or 
cons  t  ruct i on . 

The  final  report  for  this  project  is  in  three  volumes  as  follows: 

Volume  1  -  "Analytical  Modeling  of  Rock  -  Structure 
Interaction."  Final  Report,  April  1973. 

Volume  2  -  "Users  Guide  for  a  Computer  Program  for  Analytic 
Modeling  of  Rock  -  Structure  Interaction."  Final 
Report ,  April  1973  . 

Volume  3  “  "Computer  Program  for  Analytic  Modeling  of  Rock  - 
Structure  Interaction."  Final  Report,  April  1973. 

All  three  volumes  were  prepared  by  Agbabian  Associates,  under  contract 
H022CP35  with  U.  S.  Bureau  of  Mines.  The  project  was  sponsored  by  ARPA  under 
ARPA  Order  Number  1579,  Amend.  No.  3,  Program  Code  2F10. 

During  this  contract,  work  has  been  aimed  at  producing  a  user- 
oriented  computer  program.  The  work  of  writing  the  program  was  divided  into 
three  areas: 

a . 

b. 

c . 


I  nput 

Execution  and  output 
Material  properties 
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The  Input  Section  automatically  generates  the  continuum  part  of  the  finite 
element  mesh,  including  joint  el  erne,  's  ,  allows  the  user  to  add  other  elements 
(beam,  shell,  truss)  to  the  mesh,  plots  the  result,  reduces  the  bandwidth  and 
reads  loads,  material  properties,  and  other  quantities  necessary  to  the  cal¬ 
culation.  The  Execution  Section  forms  the  global  stiffness  matrix  and  selves 
equations  of  equilibrium  for  displacements  by  an  implicit  method.  The  material 
properties  are  represented  Ky  subroutines  within  the  Execution  Section,  which 
aie  written  in  a  modular  form  so  that  if  the  general  equations  of  nonlinear 
elasticity,  viscoelasticity,  viscoplasticity,  or  plasticity  do  not  suit  a 
particular  problem  they  may  be  easily  modified. 

One  of  the  guidelines  for  this  project  was  to  consolidate  existing 
finite  element  technology  into  a  single,  general  purpose  computer  program. 
Accordingly,  the  program  uses  existing  finite  elements,  a  proven  form  of 
the  equation  of  equilibrium,  existing  material  property  descriptions,  and 
an  existing  bandwidth  reducer.  However,  a  small  amount  of  new  work  was  done. 

A  new  joinc  element  was  developed  and  an  existing  concept  for  automatic  mesh 
generation  was  greatly  extended.  Also,  a  form  of  Choleski  decomposition  was 
modified  for  efficient  use  of  mu  1 1 i buf fer i ng  ,  resulting  in  substantial 
improvement  in  efficiency  of  solving  equations  of  equilibrium  using  peripheral 
storage. 

During  this  study,  some  new  work  was  reported  by  other  ARPA/Bureau 
of  Mines  contractors  which  has  been  incorporated  in  the  program.  Among  these 
are  some  creep  data  obtained  by  W.  A.  Wawersik  of  the  University  of  Utah  and 
strength/deformabi 1 i ty  data  for  faults  by  R.  E.  Goodman,  F.  E.  Heuze,  and 
Y.  Ohnishi  of  the  University  of  California,  Berkeley. 

The  technical  work  reported  in  Volume  1  is  divided  into  two  main 
parts.  The  first  part  describes  the  range  of  possible  application  of  the 
computer  program  to  mining  engineering  problems  and  reports  an  extensive 
analytic  study  of  the  Caladay  Project  hoist  room  with  which  experimental 
data  are  compared.  The  second  part  discusses  various  aspects  of  the  computer 


i  v 
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program  and  its  theoretical  basis.  Attention  is  given  to  uhe  processing  of 
input  data  and  to  options  available  to  the  user  for  mesh  generation,  sequential 
excavation  or  construction,  automatic  bandwidth  reduction  and  plotting  of  the 
mesh.  Example  problems  which  have  been  solved  during  checkout  o<:  the  program 
are  described. 

Volume  2  explains  the  overall  operation  of  the  computer  program, 
defines  how  input  should  be  submitted  to  the  program,  defines  the  material 
property  models  in  detail,  explains  the  data  stored  on  each  tape,  and  gives 
flow  charts  of  some  subroutines.  This  input  definition  is  for  the  version 
of  the  rogram  which  operates  on  the  U  N I  VAC  1108.  This  program  is  the 
iriniury  content  of  Volume  3,  and  is  available  on  seven  track  magnetic  tape 
from  DDC-TC,  U.  S.  Department  of  Commerce,  Springfield,  Virginia  22151. 
telephone  AC  (703)  321-8517.  The  tape  is  unlabeled,  even  parity,  external 
BCD  and  is  written  at  556  BP  I .  The  tape  has  constant  record  sizes,  each 
record  being  1920  characters  long  (2*1-80  column  card  images  per  record). 

An  end  of  file  mark  follows  the  last  record  on  the  tape. 

Arrangements  to  obtain  copies  of  Volumes  2  and  3  may  also  be  made 
through  Agbabian  Associates,  250  North  Nash  Street,  El  Segundo,  California 
902*15,  telephone  AC  (213)  6*10-0576. 

This  contract  was  monitored  by  Dr.  William  J.  Karwoski ,  Spokane 
Mining  Research  Center,  U.  S.  Bureau  of  Mines. 
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INTRODUCTION 

The  purpose  of  this  contract  is  to  combine  a  number  of  recent 
advances  in  finite  element  theory  and  computer  technology  into  a  computer 
program  for  analyzing  structures  and  cavities  in  rock.  This  computer  program 
applies  to  general  three-dimensional  structures,  considers  nonlinear  material 
properties  including  homogeneous  deformation  and  inhomogeneous  deformation  due 
to  joints,  anisotropic  and  time-dependent  material  properties,  gravity  loading, 
and  sequence  of  construction  or  excavation.  Since  the  program  is  intended  for 
practical  analysis  and  design,  g^eat  effort  has  been  made  to  foresee  diffi¬ 
culties  in  using  it.  For  example,  much  tedious  work,  which  formerly  was  done 
by  the  user,  has  been  eliminated  by  sophisticated  mesh  generators  and  a  band¬ 
width  reducer.  Also,  since  many  prospective  users  may  have  access  to  small 
or  medium-sized  computers  but  may  still  wish  to  solve  large  problems 
(^000-6000  equations),  the  program  uses  up-to-date  mu  1 1 i buf fe r i ng  techniques 
for  accessing  peripheral  storage  units,  thus  dramatically  reducing  computer  run 
time  for  out-of-cors  problems.  Finally,  an  attempt  is  made  to  lengthen  the 
useful  life  of  the  program  by  making  it  simple  to  add  new  elements  and  to 
expand  the  material  property  description  and  by  making  the  program  efficient 
for  and  compatible  with  a  wide  variety  of  computers. 

The  capabilities  and  limitations  of  the  computer  program  are 
summarized  as  follows: 

a.  Two-  and  three-dimensional  geometry 

b.  Small  deformations 

c.  Inhomogeneous  material  properties 

d.  Anisotropy 
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e.  Material  nonlinearity  including  ideal  and  strain-hardening 
plasticity  and  variable  moduli  models 

f.  Viscoelasticity  and  vi scop  1  as t i ci ty 

g.  Self-loading  due  to  gravity 

h.  Sequence  of  excavation  and  construction 

i .  Stat i c  live  load i ng 

The  purpose  of  this  report  is  to  discuss  what  the  computer  program 
does  and  to  describe  its  application  to  a  practical  problem  in  mining  engineer¬ 
ing.  Much  of  the  description  of  the  program  is  given  from  the  standpoint  of 
the  prospective  user.  Mesh  generation,  the  tyres  of  elements  available,  the 
types  of  loading  and  construction  which  may  be  done  and  the  properties  of 
continuous  material  and  joints  which  are  available  are  described.  The  theory 
underlying  the  present  finite  element  formulation  is  described.  The  structure 
of  the  code  is  indicated  by  logic  diagrams. 

The  application  of  the  program  to  a  practical  mining  situation  is 
reported  from  the  standpoint  of  a  user  approaching  such  a  problem  for  the 
first  time.  The  geologic  conditions  of  the  area,  and  the  joint  and  fault 
pattern  near  the  cavity  which  is  being  excavated  are  described.  The  candidate 
material  properties  are  discussed.  Calculations  are  made  using  candidate 
stress/strain  properties  and  compared  with  measurements  of  wall  deflections 
in  an  auxiliary  drift.  Parameters  defining  the  stress/strain  relations  are 
varied  until  agreement  between  calculation  and  measurement  is  reached  for 
the  auxiliary  drift.  The  same  properties  are  used  to  analyze  the  three- 
dimensional  response  of  a  cavity  in  a  similar  geologic  formation  to  excavation 
considering  creep  of  the  surrounding  rock. 
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SECTION  2 

CAPABILITY  Gi-  COMPUTER  PROGRAM  AND  ILLUSTRATIVE  CASE  HISTORY 

This  section  mentions  the  physical  aspects  of  rock  and  support 
systems  which  are  represented  by  the  computer  program  and  illustrates  part 
of  this  capability  by  analyzing  the  excavation  of  a  chamber.  The  main 
purpose  of  this  section  is  to  show  the  scope  of  problems  which  the  program 
is  capable  of  analyzing  and  to  illustrate  how  various  options  in  the  program 
are  selected  to  solve  a  specific  problem. 

The  main  work  reported  in  this  section  is  an  analysis  of  a  hoist 
room  roof  at  the  Caladay  project,  Osburn,  Idaho,  which  is  excavated  in 
argillaceous  quartzite.  Deflections  which  were  measured  during  part  of  the 
excavation  are  compared  with  the  results  of  the  analysis.  One  goal  in  per¬ 
forming  the  analysis  is  to  obtain  as  close  an  agreement  as  possible  using 
the  available  data  on  sequence  of  construction,  properties  of  rock  and  in 
situ  stresses.  The  main  benefit  of  analyses  such  as  these  performed  for  the 
Caladay  hoist  room  is  the  insight  which  they  provide  into  the  pattern  of 
stress  and  deformation.  This  insight,  which  '.he  present  computer  program 
brings  within  reach  of  even  modest  computing  facilities,  has  a  role  to  play 
a  longs i de  empi r i cal  design,  field  measurements,  and  inspection. 

2 • 1  PHYSICAL  ASPECTS  OF  ROCK  AND  SUPPORT  SYSTEMS  WHICH  THIS  PROGRAM  REPRESENTS 

To  illustrate  the  capability  of  the  program,  three  hypothetical 
problems  are  described  below.  These  problems  have  not  actually  been  solved 
with  the  present  program,  but  they  could  be  solved  at  any  time.  The  first 
is  illustrated  in  Figure  2-1.  A  section  of  tunnel  is  to  be  excavated  in  a 
region  containing  a  major  joint.  The  properties  of  the  joint  are  assumed  to  be 
known.  The  rock  adjacent  to  the  joint  is  assumed  to  be  homogeneous  and  to 
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have  viscous  properties  which  can  be  represented  by  a  visco-plastic  model. 

In  Step  1,  the  tunnel  has  not  yet  been  excavated.  Stress  in  the  rock  is 
computed  by  applying  static  overburden  to  the  edges  of  the  finite  element 
mesh.  Then  the  tunnel  is  excavated  by  removing  appropriate  elements.  At 
each  stage  of  the  excavation,  the  tunnel  roof  is  propped  by  truss  and  beam 
elements.  Eventually,  the  tunnel  is  fully  open  and  the  final  supports  are 
installed.  Each  stage  is  associated  with  an  elapsed  time,  during  which  the 
rock  flows  in  a  visco-plastic  manner.  At  each  intermediate  stage  and  at  a 
stage  the  user  defines  to  be  final,  the  stresses  in  the  rock  and  in  the  support 
elements  are  printed. 


The  second  problem  is  illustrated  in  Figure  2-2.  A  bank  is  to  be 
excavated  in  a  rock  such  as  shale  having  nonlinear,  anisotropic  stress/strain 
properties  and  an  anisotropic  fracture  criterion.  In  Step  1,  the  in  situ 
states  of  stress  are  computed  by  applying  gravitational  forces  in  a  step-by- 
step  fashion  throughout  the  grid.  In  subsequent  steps,  elements  are  removed 
in  any  sequence  the  user  desires.  Between  excavation  steps  the  remaining 
rock  will  be  checked  for  fracture  which  would  correspond  to  spalling  and 
sliding  in  an  actual  field  situation. 


The  third  example  is  an  extension  to  three  dimensions  of  the  first 
example.  The  final  stage  in  the  calculation,  at  which  the  section  of  tunnel 
under  consideration  is  fully  excavated,  is  illustrated  in  Figure  2-3.  It 
was  intended  to  develop  a  three-dimensional  joint  element  as  part  of 
the  present  project.  The  theoretical  aspects  of  this  work  are  about  50  percent 
completed  as  of  the  date  of  this  report. 


2-2  ILLUSTRATIVE  CASE  H I ST0RY--CALADAY  HOIST  ROOM 

An  analysis  was  made  of  the  excavation  of  a  hoist  room  at  the  Caladay 
Project  near  Osburn,  Idaho.  The  location  of  the  Caladay  project  is  illustrated 
in  Figure  2-A.  The  geology  of  the  area  may  be  described  as  follows: 
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2.2. 1  GEOLOGIC  DESCRIPTION 

Figure  2-4  shows  the  general  location  of  the  Caladay  underground 
hoist  room  project  in  the  Coeur  d'Alene  mining  district  in  Idaho.  The 
Coeur  d'Alene  mining  district  (Reference  2-1),  one  of  the  preeminent  lead-, 
zinc-,  and  silver-producing  areas  of  the  world,  is  near  the  base  of  the 
panhandle  of  northern  Idaho.  Spokane,  Washington,  is  about  75  miles  to  the 
west,  and  Missoula,  Montana,  about  MO  miles  to  the  east  of  the  district. 

The  district  lies  wholly  within  the  Coeur  d'Alene  Mountains  which  are  a 
rugged,  deeply  dissected  mountain  mass.  Rock  units  of  the  Coeur  d'Alene 
district  are  metamorphosed  sedimentary  of  the  Precambrian  Belt  Series.  This 
series  consists  of  six  formations.  The  three  of  interest  are  the  Revett 
Quartzite,  St.  Regis,  and  Wallace  Formations.  They  are  in  a  complex  anti- 
cl ine-syncl ine  structural  system.  A  very  complex  geologic  history  associated 
with  strike  slip  movements  along  the  Osburn  fault  system  has  caused  individual 
beds  to  be  sharply  contorted,  crumpled,  overturned,  and  faulted.  Formation 
contacts  are  not  well  defined. 

The  Revett  Quartzite,  St.  Regis,  and  Wallace  are  three  of  the  oldest 
formations.  The  St.  Regis  forms  its  contact  with  the  Revett  Quartzite  Forma¬ 
tion  and  grades  from  interbedded  quartzite  and  argillite  upward  into  a  dominantely 
thin-bedded  argillite.  The  upper  several  hundred  feet  of  the  St.  Regis  are 
characterized  by  a  finely  laminated  argillite.  The  Wallace  Formation  is  a 
heterogeneous  group  of  rock  comprising  quartzite,  argillite,  dolomite,  and 
1 imestone. 

The  rock  of  the  Belt  series  is  fine  grained,  and  the  principal 
minerals  are  quartz  and  sericite  with  quartz  the  more  abundant  mineral. 

Faults  are  the  dominant  structural  features  (Figure  2-4)  in  the 
area,  and  a  myriad  of  them  cut  the  country  rock  into  a  complex  pattern  of 
blocks.  Hundreds  can  be  seen  underground.  As  a  result  of  the  regional 
faulting  and  folding,  most  of  the  beds  dip  greater  than  45  deg.  The  rock 
of  the  district  was  intensely  deformed  in  a  complex  pattern  which  can  he 


8 


FIGURE  2-4.  FAULT  LOCATIONS,  COEUR  D'ALENE  DISTRICT 


R-721 5-1-2701 


R- 7215-1-2701 


referred  to  as  a  structural  knot.  Much  of  this  regional  geology  is  similar 
to  the  local  geology  at  the  Caladay  project  site. 

The  Caladay  project  (Reference  2-2)  was  developed  by  driving  an  adit 
5,200  ft  into  a  mountainside.  The  adit  begins  in  the  Wallace  Formation.  At 
2,900  ft  the  Wallace/St.  Regis  contact  is  found.  The  adit  intersects  the 
Polaris  fault  at  approximately  3,200  ft  where  the  Wallace  Formation  is 
reentered.  The  Wallace/St.  Regis  contact  is  encountered  again  at  A, 800  ft. 

The  hoist  room  is  about  A00  ft  from  this  Wallace/St.  Regis  contact.  Thus/ 
the  hoist  room  is  in  the  upper  several  hundred  feet  of  the  St.  Regis  Formation. 

Bedding  in  the  hoist  room  dips  nominally  65  deg  and  strikes  normal 
to  the  long  axis  of  the  room.  Bedding  in  the  hoist  room  can  be  classed  into 
thee  broad  categories  based  on  its  mineral  composition.  Zone  A  is  estimated 
to  contain  60  percent  argillite  and  A0  percent  argillaceous  quartzite.  It 
is  thin  bedded  with  beds  varying  from  ]~ss  than  1  in.  up  to  6  in.  in  thickness. 
Zone  B  is  estimated  to  contain  50  percent  argillaceous  quartzite,  A 0  percent 
argillite,  and  10  percent  quartz.  The  quartz  content  diminishes  from  bottom 
to  top  of  the  room.  The  bedding  i s  A  to  8  in.  in  thickness.  Zone  C  is  prin¬ 
cipally  argillite  with  minor  amounts  of  argillaceous  quartzite,  and  this  is 
generally  thinly  bedded. 


Most  of  the  faults  in  the  area  are  slips  along  bedding  planes. 
These  faults  contain  a  silt  or  silty  clay-like  gouge.  One  predominant  fault 
(sheave  room  fault)  projects  over  the  top  corner  of  the  hoist  room  and  into 
the  sheave  room.  This  fault  is  about  2  to  10  in.  in  thickness  and  does  not 
appear  to  have  asperities  of  any  consequence.  At  the  other  end  of  the  hoist 
room  exists  a  large  cross-bedding  fault.  This  fault  dips  85  deg  across  the 
bedding  at  85  deg.  Its  shear  zone  is  estimated  at  no  less  than  50  ft  thick. 

The  gouge  material  ranges  from  large  chunks  of  broken  rock  to  sand  and  silt¬ 
like  materials. 
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2.2.2  PHYSICAL  SITUATION  AND  MATHEMATICAL  MODEL  FOR  CALADAY  HOIST  ROOM 

Figure  2-5  shows  the  drift  and  tunnel  complex  surrounding  the  hoist 
room  at  the  Caladay  project.  The  figure  also  shows  the  local  geologic  struc¬ 
ture  including  distribution  of  different  kinds  of  rocks,  location,  and 
orientation  of  bedding  faults  and  cross-cutting  major  faults.  The  two  major 
types  of  rock  encountered  in  the  area  are  argillite  and  quartzite.  These  are 
found  mostly  as  argillaceous  quartzite.  A  major  cross-bed  fault,  located 
next  to  the  southwest  wall  of  the  hoist  room,  contains  soft  clay  and  shattered 
argillite. 

Some  details  concerning  the  hoist  room  are  given  in  Figure  2-6. 

The  excavation  of  the  room  was  started  from  the  top  of  the  room,  which  was 
reached  through  a  raise  dug  upward  from  the  entry  drift.  The  excavation  was 
then  carried  downward  by  presplitting  and  mucking  horizontal  layers  of  material 
until  the  bottom  of  the  room  was  reached. 

To  monitor  the  response  of  the  back  of  the  hoist  room,  three  exten 
someters  of  30- f t  length  were  installed  in  the  back  as  soon  as  it  was  reached 
at  the  top  of  the  access  raise.  Readings  on  these  meters  in  the  course  of 
excavation  were  subsequently  taken  for  a  period  of  approximately  seventy  days 
until  access  to  these  meters  was  lost  at  the  end  of  the  top  slice  and  first 
bench  excavation.  These  readings  were  later  processed  into  the  form  of  strain 
rates  and  cumulative  displacements  between  anchors  located  along  the  length 
of  the  extensometers .  Figure  2"7  shows  the  time  variation  of  cumulative  dis 
placement  between  the  toe  and  collar  of  each  one  of  the  extensometers.  One 
of  the  objectives  of  the  finite  element  analyses  to  be  presented  herein  is 
to  check  on  these  observed  values. 

To  introduce  the  modeling  procedure,  one  of  the  preliminary  two- 
dimensional  finite  element  meshes  which  is  used  only  to  investigate  mesh 
fineness,  loading  and  boundary  supporting  conditions  suitable  for  eventual, 
three-dimensional  modeling  of  the  hoist  room  is  shown  in  Figure  2-8(a).  This 
11  by  11  mesh  was  automatically  generated  and  plotted  by  the  code  from  the 
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input  key  diagram  shown  in  Figure  2-8(b).  The  so-called  "key  diagram"  is 
discussed  in  detail  in  Section  5.  Figure  2-8(c)  shows  the  configuration  of 
the  key  diagram  when  the  actual  coordinates  of  the  controlling  nodal  points  are 
taken  into  account.  The  key  diagram  has  been  so  planned  that  Elements  39,  SO, 
61,  72,  83,  and  Elements  40,  51,  62,  73,  and  84  outline  the  top  slice  and  the 
first  bench,  respectively,  and  can  be  removed  in  sequence  in  the  course  of 
numerical  calculation  to  simulate  the  actual  excavation  sequence. 

The  mathematical  model  is  roller-supported  along  the  left-  and 
right-hand  edges  except  at  Nodal  Points  7  and  139  which  are  hinged.  A 
uniform  unit  pressure  (l  ps i )  is  applied  along  the  top  and  bottom  edges  to 
simulate  the  overburden  pressure.  Since  the  model  is  elastic,  the  response 
to  any  other  vertical  pressure  can  be  readily  computed  by  scaling.  These 
loading  and  supporting  conditions  are  the  result  of  earlier  comparison  studies 
which  show  that,  for  this  plane  strain  idealization,  applying  the  overburden 
pressure  along  two  edges  (top  and  bottom)  is  superior  to  applying  it  along 
one  edge  only  (say,  the  top  edge)  and  supporting  the  other  edge  (the  bottom 
edge)  as  well.  However,  applying  normal  pressure  to  each  of  three  orthogonal 
faces  to  represent  in  situ  stresses  is  more  suitable  for  the  refined  three- 
dimensional  model  described  below. 

Figure  2-9  summarize  the  distribution  of  vertical  normal  stress 
determined  from  the  preliminary  plane  strain  analysis  with  idealized  elastic 
properties.  The  narrow  layer  of  argillaceous  quartzite  is  assumed  to  be 
approximately  twice  as  stiff  as  the  argillaceous  quartzite  on  the  left  and 
twenty  times  as  stiff  as  the  argillite  on  the  right.  Figure  2-9  (a)  shows 
the  stress  distribution  prior  to  excavation.  Figures  2-9(b)  and  2-9 (c)  show 
the  distribution  following  the  excavation  of  the  top  slice  and  the  first 
bench,  respectively.  The  latter  solutions  are  the  result  of  performing  four 
iterations  following  the  reformulation  of  the  global  stiffness  matrix  to 
delete  the  stiffness  of  the  elements  corresponding  to  the  part  of  the  room 
newly  excavated.  Although  comparison  with  a  closed  form  solution  cannot 
be  made  ror  the  case  of  inhomogeneous  rock,  the  familiar  pattern  of  stress 
concentration  around  a  rectangular  hole  is  clearly  seen.  No  attempt  is 
made  to  compare  the  results  of  this  computation  with  the  extensometer 
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(a)  DISTRIBUTION  OF  VERTICAL  NORMAL  STRESS  BEFORE  EXCAVATION 
(ALL  STRESSES  ARE  COMPRESSIVE) 


(b)  DISTRIBUTION  OF  VERTICAL  NORMAL  STRESS  FOLLOWING  TOP-SLICE  EXCAVATION 
(ALL  STRESSES  EXCEPT  THOSE  WITH  POSITIVE  SIGNS  ARE  COMPRESSIVE) 
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DISTRIBUTION  OF  VERTICAL  NORMAL  STRESS  FOLLOWING  FIRST-BENCH  EXCAVATION 
(ALL  STRESSES  EXCEPT  THOSE  WITH  POSITIVE  SIGNS  ARE  COMPRESSIVE) 


FIGURE  2-9.  RESULTS  OF  PRELIMINARY  TWO-DIMENSIONAL  ELASTIC  ANALYSIS 


R-721 5- 1 “2701 

data  because  (a)  the  material  coefficients  used  in  this  early 
study  differ  from  those  finally  adopted,  (b)  the  major  fault 
described  earlier,  which  would  appear  off  the  left  end  of  the  hole  in  the 
current,  two-d imens ional  mesh,  is  not  taken  into  account  In  this  model,  and 
(c)  a  two-dimensional,  plane  strain  analysis  in  this  case  is  not  expected  to 
yield  a  good  solution  especially  when  the  plane  chosen  is  parallel  to  the 
long  axis  of  the  hoist  room.  The  solution  is  presented  here  as  an  introduction 
to  the  more  complicated  analyses  which  follow. 

2.2.3  MATERIAL  PROPERTIES  FOR  ANALYSIS  OF  CALADAY  PROJECT  HOIST  ROOM 

This  section  describes  properties  of  the  St.  Regis  formation 
argillaceous  quartzite  and  argillite  which  are  used  in  computing  stresses 
and  deformations  at  the  Caladay  hoist  room.  The  goal  of  defining  the  properties 
is  to  determine  the  elastic  and  inelastic  stiffnesses  and  strengths  which 
are  typical  of  rock  masses  In  situ.  In  the  analysis  it  is  assumed  that  these 
properties  are  homogeneous  throughout  the  volume  of  a  finite  element,  which 
in  the  present  case  is  about  10  to  15  feet  on  a  side. 

The  material  properties  are  selected  from  laboratory  and  in-situ 
measurments  and  by  consultation  with  geologists  and  mining  engineers  who  are 
familiar  with  the  site.  Due  to  wide  variations  in  the  types  of  rock  present 
and  to  uncertainty  in  translating  laboratory  properties  into  in-situ  properties, 
corsultation  is  considered  to  be  the  most  important  element  in  selecting 
reasonable  properties. 

The  fundamental  decision  which  was  made  regarding  material  properties 
is  the  mathematical  idealization  within  which  the  properties  are  represented. 

As  explained  in  subsequent  sections  of  this  report,  the  available  models 
include  such  features  as  anisotropy,  plasticity,  viscoelasticity,  and 
viscoplasticity.  The  decision  to  use  isotropic  viscoelasticity  is  based 
partly  on  positive  factors,  such  as  recommendations  from  knowledgeable 
geologists  and  laboratory  and  in  situ  measurements,  and  partly  on  negative 
factors  such  as  the  absence  of  sufficient  data  on  which  to  base  an  anisotropic 
model.  Viscoplasticity  is  an  attractive  model  for  argillite  and  argillaceous 
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quartzite  because  both  plasticity  and  time  dependence  appear  to  be  present. 
Precedent  for  using  plasticity  to  represent  these  materials  is  found  in 
Reference  2-3.  Clear  evidence  of  time-dependent  deformations  is  also  found 
in  this  and  in  similar  work  (Reference  2-4).  The  viscoplastic  model  was 
rejected  primarily  for  lack  of  data  which  is  required  to  distinguish  the 
amount  of  t i me- i ndependen t  inelasticity  from  time-dependent  inelasticity. 
Without  rational  basis  for  such  a  distinction  it  is  better  to  use  a  simpler 
model  with  only  one  mechanism  for  inelasticity.  Although  an  elastic/plastic 
model  is  used  for  the  analyses  of  Reference  2-3,  the  present  case  involves 
lower  in  situ  stresses  and  hence  a  smaller  contribution  to  displacements  from 
time-independent  plasticity.  Based  on  this  reasoning,  time-dependent  elastic 
properties,  as  are  represented  by  linear  viscoelasticity,  are  assumed  to  be 
the  dominant  mechanism  governing  stresses  and  deformations  around  the  hoist 
room.  The  model  of  stress/strain  properties  for  argillite  is  developed 
according  to  this  assumption. 

In  contrast  to  the  argillite  and  argillaceous  quartzite,  the  cross 
bed  fault  which  passes  near  the  hoist  room  is  filled  with  clay  minerals  and 
fractured  argillite.  The  shear  strength  of  this  material  has  not  been 
measured  due  to  the  difficulty  of  sampling  such  an  inhomogeneous  material. 

It  is  likely,  however,  that  the  shear  strength  of  wet  fault  gouge  is  very 
low.  Accordingly,  it  would  be  well  represented  by  an  elastic/perfectly  plastic 
model.  The  principal  numerical  difficulty  in  using  such  a  model  arises  from 
applying  stress  boundary  conditions,  which  in  the  present  analysis  represent 
in  situ  stresses.  A  simple  calculation  shows  why  it  is  not  possible  to 
specify  both 

a.  Shear  strength  of  fault  gouge  of  order  100  psi  and 

In  situ  stress  components  of  1200,  3000,  and  3500  psi 

The  maximum  shear  stress  associated  with  the  in  situ  stresses  is  (3500  psi  - 
1200  psi)/2  =  1150  psi,  which  is  much  larger  than  any  reasonable  shear 
strength  for  the  fault  gouge.  Clearly  the  values  of  in  situ  stress  which  are 
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shown  in  Figure  2-6  are  average  or  representative  values  which  do  not  apply 
near  the  fault.  The  following  are  three  alternatives  in  representing  the 
fault  within  the  capability  of  the  present  computer  program. 

a.  Move  the  boundaries  of  the  mesh  far  beyond  the  fault,  use 
elastic/plastic  properties  for  the  fault  ann  allow 
plastic  flow  to  modify  the  local  stress  field  as  occurs 
in  nature. 

b.  Use  very  low  elastic  moduli  for  the  fault  material,  thus 
allowing  the  stiffer  argillite  and  argillaceous  quartzite 
to  sense  a  flexible  zone  without  interfering  with  the 
stress  boundary  conditions. 

c .  Omit  the  fau  1 1 . 

Alternative  (a)  is  rejected  on  account  of  cost.  Alternative  (b)  is 
superior  to  Alternative  (c) ,  aid  hence  is  selected. 

2.3.3. 1  Viscoelastic  Properties  of  Argillite 

The  viscoelastic  model  for  argillite  (zone  c)  is  based  on  two 
assumptions  which  are: 

a.  The  argillite  may  be  represented  by  three  parameter  solid 
containing  a  Kelvin  element  (spring  and  dashpot  in  parallel) 
in  series  with  a  spring. 

b.  The  elastic  bulk  and  shear  moduli  of  the  in  situ  argillite 
rock  mass  can  be  found  on  the  basis  of  laboratory  tests 
and  judgment. 

The  viscoelastic  properties  are  selected  such  that  analysis  of  deformation 
around  an  auxiliary  shaft  matches  the  field  data.  In  this  way  the  model 
incorporates  field  and  laboratory  measurements  and  engineering  judgment. 
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TABLE  2-1.  REPRESENTATIVE  ELASTIC  PROPERTIES  OF  PURE  MINERALS 
BASED  ON  LABORATORY  MEASUREMENTS 


Type 


Pure  Quar tz i te 
Pure  Arg i 1 1 i te 


In  Situ  Young's  Modulus 

In  Situ  Poisson's  Ratio 

7.0  x  106  psi 

0.2 

0.35  x  106  psi 

0.1 

converted  to  K  and  G  as  follows: 

K  -  TTi^T 

(2-3) 

°  ■  2li^ 

(2-i*) 

The  viscoelastic  coefficients  to  be  determined  for  substitution  into 
Equations  2-1  and  2-2  are  l/nK.  1/nG  and  the  K  and  G  for  the  Kelvin 
element.  This  is  done  by  parametric  variation  of  these  coefficients  in  a 
finite  element  analysis  to  match  an  in  situ  experiment. 

The  experiment  is  illustrated  in  Figure  2-10  and  is  described  in 
Reference  2-3.  At  a  typical  station,  rock  bolts  were  installed  as  close  to 
the  advancing  face  as  possible.  The  distances  BB1  and  CC 1  were  measured.  As 
excavation  of  the  drift  progressed,  movements  of  the  rock  bolts  were  monitored 
by  means  of  extensome ters .  The  results  of  these  measurements,  reproduced  from 
Reference  2-3,  are  shown  in  Figure  2-11.  The  viscoelastic  coefficients  were 
selected  by  assuming  values  for  and  Dq ,  applying  the  measured  values 

of  in  situ  stresses  shown  in  Figure  2-12  to  the  finite  element  mesh  in  the 
same  figure,  then  changing  the  values  until  the  agreement  between  measurement 
and  calculation  is  satisfactory.  In  practice,  the  argillite  :n  Zone  C  was 
first  assumed  to  be  elastic  in  bulk  and  viscoelastic  in  shear.  Then  it  was 
assumed  to  be  viscoelastic  in  bulk  and  elastic  in  shear.  It  was  originally 
intended  to  assume  combined  viscoelasticity  in  shear  and  bulk.  However,  visco¬ 
elasticity  in  bulk  alone  seemed  to  be  adequate.  The  viscoelastic  coefficients 
finally  selected  are  shown  in  Table  2-2.  The  comparison  between  in  situ 
measurements  and  the  finite  element  calcualtions  using  these  properties  is 
shown  in  Figure  2-11. 

23 


— — - — - - 


~  . . .  — 


R-721 5-1-2701 


(a)  LOCATION  OF  ROCK  BOLTS  SERVING  AS  MEASURING  POINTS  FOR  MEASURING 
RELATIVE  DISPLACEMENTS  IN  DRIFT  USING  EXTENSOMETER 


FIGURE  2-10.  AUXILIARY  DRIFT  IN  SILVER  SUMMIT  MINE  USED  TO  DETERMINE 
IN  SITU  VISCOELASTIC  PROPERTIES  OF  ARGILLITE 
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CREEP  MODEL  FOR  ARGILLITE 


/ 


(b)  A  REPRESENTATIVE  CROSS-SECTIONAL  VIEW  AT  STATION  SS- 3  OF  THE 
SOUTH  EAST  LATERAL,  AOOO  LEVEL,  SILVER  SUMMIT  MINE 
(REFERENCE  2-5) 


FIGURE  2-10.  (CONTINUED) 
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The  stress/strain  relations  for  the  Kelvin  portion  of  the  visco¬ 
elastic  model  are  divided  into  two  parts.  One  part  governs  volumetric 
behavior  and  a  second  part  governs  shear  behavior  as  follows: 


0 . 


1  1 


— 5  =  K  e .  .  +  n 


e .  . 


K  1  i 


=  2G  e!.  + 
U 


nr  e : . 

G  1  j 


(2-1) 

(2-2) 


where 


e .  . 

i  i 

“ 

Volumetric  strain 

1 

3  i  i 

= 

Mean  normal  stress 

e !  . 

'J 

= 

Deviatoric  strain  component 

o!  . 

'J 

= 

Deviatoric  stress  component 

k,g 

= 

Elastic  bulk  and  shear  moduli 

VnG 

= 

Viscous  coefficients  in  bulk  and 

and  nG 

are 

independent  coefficients.  Thus, 

both  components. 


The  first  step  in  developing  the  model  for  the  aroillite  is  to 
determine  in  situ  elastic  moduli.  Laboratory  measurements  of  Young's  moduli 
and  Poisson's  ratios  are  available  for  argillites,  quartzites  and  mixtures 
of  the  two  in  Reference  2-6.  Host  of  these  data  are  for  rocks  from  the  St. 
Regis  and  Revett  Quartzite  formations.  The  lowest  modulus  for  pure  quartzite 
■s  about  4  x  10  psi  and  the  highest  is  about  13  x  106  psi.  Poisson’s  ratio 
appears  to  vary  from  about  0.05  to  about  0.24.  Based  on  these  measurements 
and  suggestions  received  in  Reference  2-7,  the  properties  for  pure  quartzite 
are  assumed  to  be  as  shown  in  Table  7-1.  Based  on  suggestions  in  Reference  2-7 
the  in  situ  moduli  of  the  argillaceous  quartzites  in  the  vicinity  of  the  hoist 
room  are  computed  from  the  fractions  of  pure  argillite  and  pure  quartzite. 

These  are  indicated  in  Table  2-2. 


29 


aiaaaiiMiaii 


I  iriiihiatii 


Jl  .  -I  ,1.  II  llH^pp 


i^i^liiiWii^H  UWWWjP^y^fiiyajWjy  li  J.^.,1..  IIIL1  j  jljp  .UniWg.  I  .. 


A 


,i  .mu 


R-721 5- 1-2701 


2. 2. 3. 2  Elastic  Properties  of  Argillaceous  Quartzite  and  Fault  Zone 

The  properties  of  the  argillaceous  quartzite  in  Zones  A  and  B  are 
assumed  to  be  linearly  elastic  as  shown  in  Table  2-2.  The  elastic  moduli  are 
determined  from  the  percentage  of  quartzite  and  argillaceous  quartzite  in 
each  as  described  above.  The  judgment  of  experienced  engineers  who  are 
familiar  with  the  area  is  the  most  important  factor  in  selecting  these 
coefficients.  It  would  have  been  preferable  to  include  viscoelasticity  in 
the  models  of  material  in  Zones  A  and  B.  However,  it  was  not  possible  to 
obtain  suitable  samples  of  the  rock  in  time  for  the  present  analysis.  More¬ 
over,  measurements  on  rock  which  was  similar  but  which  contained  more  quartzite 
showed  no  appreciable  creep. 

The  moduli  of  material  in  the  fault  was  arbitrarily  assumed  to  be 
very  low,  as  shown  in  Table  2-2. 


2.  3  THREE-DIMENSIONAL  ANALYSIS  OF  HOIST  ROOM  AT  CALADAY  PROJECT 

The  physical  situation  of  the  hoist  room  which  is  shown  in  Figure  2-6 
and  the  mathematical  models  for  the  various  types  of  rocks  which  are  described 
above  arc  combined  into  a  three-dimensional  finite  element  analysis. 

The  finite  element  mesh  for  the  final  three-dimensional  analysis  of 
the  hoist  room  is  shown  in  Figure  2-13.  It  is  an  11  x  11  x6  mesh  with  1008 
nodal  points  and  726  elements.  It  is  loaded  on  three  of  its  six  faces  to 
simulate  the  in  situ  stress  condition,  and  simply  supported  on  the  remaining 
three  faces  as  indicated.  Except  for  points  on  the  simply  supported  faces, 
each  point  has  three  translational  degrees  of  freedom  so  that  the  simultaneous 
solution  of  2712  equation  is  involved  in  each  load  step.  The  four  idealized 
zones  having  different  materials  are  shaded  in  the  figure.  These  are  the 
same  zones  as  shown  in  Figure  2-6  and  Table  2-2. 

For  the  convenience  of  creating  plane  interfaces  among  the  different 
material  zones,  all  the  mesh  lines  running  in  the  z-direction  have  been  made 
straight  and  parallel  to  the  z-axis.  The  mesh  has  been  created  in  such  a 
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FIGURE  2-13.  THREE-DIMENSIONAL  FINITE  ELEMENT  MESH  WITH  ORIGINAL 
CONFIGURATION  OF  TOP  SLICE  AND  FIRST  BENCH  FOLLOWED 
BY  THREE  STAGES  OF  EXCAVATION 
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manner  that  the  rectangular  cavity  to  be  created  by  the  excavation  of  the  top 
slice  and  the  first  bench  is  represented  by  two  layers  of  fifteen  elements 
each.  These  thirty  elements  are  "turned  off"  in  the  course  of  the  numerical 
calculation  in  the  sequence  illustrated  in  the  detailed  drawing  to  simulate 
the  actual  excavation  sequence  that  took  place  in  the  field.  The  mesh  is 
made  finer  around  the  cavity  for  better  stress  definition  and  the  first  three 
layers  of  mesh  lines  over  the  back  are  located  at  6  ft,  15  ft,  and  30  ft  above 
the  back  to  facilitate  the  checking  of  extensometer  readings.  The  mesh  as 
designed  is  rather  coarse;  but  it  is  as  fine  as  is  permitted  by  the  computer 
budget  of  the  present  project. 

Computation  begins  with  the  mesh  representing  solid  rock  except 
for  a  small  cavity  representing  the  subdrift  which  was  driven  before  the 
extensomete rs  were  installed.  The  first  step  of  the  computation  ^involves 
applying  the  in  situ  stresses  to  the  mesh  and  finding  the  corresponding  stress 
and  displacement  fields.  The  second  step  is  an  iteration  o  obtain  an  improved 
approximation  to  equilibrium  conditions.  In  the  computation,  the  end  of 
Step  2  corresponds  to  installation  of  the  extensometers ,  and  hence  to  a  zero 
value  of  relative  displacement.  Between  Steps  2  and  3,  time  is  assumed  to 
elapse  and  viscoelastic  flow  occurs  in  the  argillite.  Subsequently,  elements 
are  turned  off  at  appropriate  intervals  to  simulate  the  progress  of  the 
excavation,  and  the  real  time  is  incremented  intermittently  to  allow  creep 
to  take  place  in  between  the  excavation  steps.  The  entire  sequence  of 
operations  is  summarized  in  Table  2-3.  The  solution  is  iterated  once  follow¬ 
ing  each  excavation  step  which  involves  the  reformulation  of  the  global  stiff¬ 
ness  matrix.  Additional  iterations  would  be  highly  desirable,  but  were 
omitted  to  keep  within  the  computer  budget. 

The  result  of  the  run  is  represented  by  the  displacement  components 
of  all  the  node  points  and  the  stress  and  strain  components  in  all  the  elements 
for  each  one  of  the  twelve  steps  taken.  It  i s  of  interest  to  examine  the 
stress  distribution  in  the  mesh  at  the  end  of  the  second  step.  Since  the 
hole  that  exists  in  the  mesh  at  this  time  is  still  small,  the  stress  distribution 
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TABLt  2-3.  TWELVE  STEPS  TAKEN  IN  THE  FINAL 
THREE-DIMENSIONAL  CALCULATION 


Rea  I  Time 
(Days) 


0 

(June  12) 
0 
13 

(July  1) 


Purpose 

To  create  the  in  situ  stress  condition  with  the 
subdrift  already  in  place. 

Iterate  to  stabilize  the  elastic  solution. 

To  allow  creep  to  take  place  for  18  days. 


Advance  top  slice  excavation  by  removing  six 
additional  elements. 


18 


To  perform  an  iteration. 


33 

(July  15) 


To  allow  15  additional  days  of  creep  to  take 
place. 


33 

33 


Remove  six  more  elements  to  complete  the  excava¬ 
tion  of  the  top  slice. 

To  perform  an  iteration. 


^9 

(August  1) 
^9 
^9 


To  allow  16  additional  days  of  creep  to  take 
p 1  ace . 

Excavate  the  entire  first  bench  in  one  step. 
To  perform  an  iteration. 

To  allow  the  final  22  days  of  creep  to  take 
p 1  ace . 


71 

(August  22) 


R-721 5-1-2701 


A 

created  in  the  mesh  should  be  very  nearly  the  in  situ  stress  condition  that 
can  be  realized  in  the  mathematical  model  under  the  loading  and  boundary 
conditions  adopted.  Figure  2-14  summarizes  the  distribution  of  the  three 
normal  stresses  in  the  direction  of  the  coordinate  axes  in  the  first  and  the 
last  layer  of  elements  when  the  mesh  is  viewed  from  the  negative  z-direction. 

It  is  noted,  first  of  all,  that  the  stress  distribution  is  very 
irregular,  with  the  stress  values  in  some  portions  of  the  mesh  quite  unlike 
the  homogeneous  stress  condition  which  is  applied  to  three  faces  of  the  mesh. 
What  is  mainly  responsible  for  this  dramatic  variation  in  the  stresses  is 
the  difference  among  the  stiffnesses  of  the  materials  involved  in  the  four 
zones.  Young's  Modulus  varies  as  much  as  two  orders  of  magnitude  from  the 
softest  to  the  stiffest  material.  Under  these  circumstances,  there  are  large 
amounts  of  stress  transfer  from  the  softer  to  the  stiffer  zones,  causing 
large  shearing  stresses  at  the  material  interfaces  and  very  uneven  reaction 
distribution  over  the  three  simply  supported  faces. 

It  is  assumed  in  the  calculation  that  the  extensometers  are  installed 
at  the  end  of  the  second  step  so  that  they  are  stretched  only  by  the  additional 
deformation  created  henceforth.  On  the  basis  of  this  assumption,  the  equival¬ 
ent  of  two  extensometer  records  have  been  deduced  from  the  computer  run  and 
are  plotted  in  Figure  2-15.  Figure  2- 1 6  shows  the  variation  of  the  vertical 
normal  stresses  over  the  back  of  the  room  as  excavation  progresses  sidewise 
from  the  subdrift  and  then  downward  into  the  first  bench.  Figure  2-1/  shows 
the  time  history  of  the  downward  deflection  of  the  center  of  the  back  end  the 
upheaval  of  the  center  of  the  floor  of  the  first  bench  during  the  same  period. 
In  calculating  the  deflections  oF  the  back  and  the  floor  during  the  excava¬ 
tions,  the  deflected  configuration  of  the  mesh  at  the  end  of  the  second  step 
is  considered  as  the  datum.  Based  on  this  datum,  the  upheaval  of  the  entire 
first  bench  floor  at  the  end  of  71  days  (12th  step)  is  pictured 
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(c)  AT  TIME  -  33  DAYS  (END  OF  STEP  8) 


ABSOLUTE  DEFLECTION  OF  CENTER  OF  BACK  AND  FIRST  BENCH  FLOOR 
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in  Figure  2- 1 8 (a ) .  Figure  2- 1 8 ( b )  shows  the  incremental  deflection  of  the 
same  floor  from  49  to  71  days  (llth  to  12th  step)  due  entirely  to  creep. 

Figure  2-l8(a)  is  of  !r,terest  only  from  an  applied  mechanics  point  of  view, 
because  it  shows  the  total  distortion  of  the  plane  corresponding  to  the  first 
bench  floor  due  to  excavation  and  creep,  including  a  part  which  occurred  before 
the  floor  is  exposed.  Figure  2- 1 8 ( b)  is  of  interest  to  the  designer  who  is 
concerned  primarily  with  the  floor  heave  which  occurs  after  the  floor  has  been 
cleared.  If,  for  example,  a  concrete  foundation  were  to  be  laid  on  the  first 
bench  floor  at  the  end  of  49  days,  its  shape  at  the  end  of  71  days  would  be 
similar  to  that  in  Figure  2- 1 8 ( b ) .  This  distortion  may  be  too  great  for 
proper  operation  of  the  hoists.  One  alternative  would  be  to  wait  until  the 
creep  rate  has  become  negligible  before  laying  the  foundation.  A  second 
alternative  would  be  to  anticipate  the  eventual  shape  and  to  compensate  by 
digging  the  floor  initially  too  deep  in  parts  where  the  floor  is  expected  to 
heave,  and  to  lay  the  foundation  when  the  floor  reaches  a  stable,  level 
condi tion. 


42 


R-721 5-1 -2701 


A 

2.4  DISCUSSION  OF  RESULTS 

There  are  three  aspects  of  the  analysis  of  the  Caladay  Project  hoist 
room  which  would  be  modified  if  a  second  calculation  were  made.  These  are 

a.  Improving  the  method  of  introducing  in  situ  stresses 

b.  Limiting  the  region  which  undergoes  time-dependent  deformation 
to  the  zone  of  influence  of  the  hoist  room 

c.  Refining  the  finite  element  mesh  in  the  vicinity  of  the  hoist 
room 

The  present  method  of  introducing  in  situ  stresses  is  illustrated 
in  Figure  2-13.  The  result  is  shown  in  Figure  2-14  and  2- 1 6  in  terms  of  spatial 
distribution  of  stresses  at  a  selected  stage  of  the  calculation.  This  method 
provides  that  the  stresses  are  equal  to  the  assumed  in  situ  stresses  only  near 
the  faces  where  they  are  applied  as  external  forces.  Throughout  the  remainder 
of  the  finite  element  model,  stresses  are  greatly  affected  by  the  relative 
stiffnesses  and  Poisson's  ratios  of  the  four  layers  of  rock.  This  nonuniformity, 
which  probably  also  occurs  in  nature  but  to  a  lesser  degree,  produces  stress 
concentrations  which  are  1/3  to  1/2  as  great  as  those  due  to  the  hole.  This 
suggests  that  by  applying  stresses  to  three  of  the  six  faces  of  the  finite 
element  model,  the  user  does  not  obtain  the  initial  stress  conditions  he 
desires.  For  highly  inhomogeneous  models,  it  would  be  better  to  introduce 
the  desired  internal  stresses  directly,  rather  than  attempt  to  find  external 
boundary  conditions  which  produce  the  desired  internal  stress  state.  A 
slight  modification  of  the  computer  program  is  needed  to  prestress  continuum 
elements . 

There  is  one  definite  weakness  in  the  present  viscoelastic  model  and 
one  possible  weakness.  The  definite  weakness  is  that  all  of  the  argillite  in 
the  finite  element  model  undergoes  viscoelastic  volume  strain.  Thus,  even  a 
finite  element  of  argi 1 1 i te  which  is  so  distant  from  the  hoist  room  that  its 
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stress  is  unaffected  by  the  excavation  nevertheless  undergoes  creep.  This 
creep  has  no  physical  counterpart  and  hence  should  not  be  present  in  the 
mathematical  model.  In  contrast,  creep  which  has  a  physical  counterpart  occurs 
in  the  vicinity  of  the  excavation.  Physically,  this  creep  is  related  to  move¬ 
ment  along  bedding  planes  which  accompanies  changes  in  the  local  stress 
condition  due  to  excavation.  This  picture  suggests  that  the  physical  situation 
would  be  more  nearly  represented  by  assigning  creep  properties  to  a  small  zone 
of  argillite  in  the  immediate  vicinity  of  the  hoist  room.  However,  a  technique 
has  not  yet  been  developed  for  manipulating  the  finite  element  computer  program 
such  that  creep  occurs  in  this  zone  only  in  response  to  changes  in  stress  from 
the  in  situ  condition  to  a  new  condition  related  to  a  gradually  enlarging 
cavi ty . 

The  possible  weakness  in  the  present  viscoelastic  model  is  that  it 
predicts  progressive  shortening  of  the  extensometers  as  time  elapses.  This 
is  visible  in  the  falling  portion  of  the  calculated  curves  in  Figure  2-15. 

It  is  not  clear  from  the  measurements  whether  there  is  time-dependent  extension 
or  contraction  of  the  extensometers.  However,  it  appears  to  the  present 
authors  that  the  measurements  indicate  time-dependent  extension.  To  change 
the  viscoelastic  model  so  that  it  agrees  with  the  apparent  creep  measurements 
requires  replacing  the  assumption  of  viscoelastic  volume  strain.  This  is 
because  the  mean  normal  stress  in  the  rock  between  the  ends  of  the  extensometer 
is  compressive  and  the  viscoelastic  volume  strain  is  also  required  to  be  com¬ 
pressive.  It  appears  possible  that  a  different  viscoelastic  mode)  can  be 
found  to  improve  agreement.  For  example,  a  Kelvin  model  in  shear  can  probably 
be  found  which  predicts  extension  of  the  extensometers.  However,  as  is  pointed 
out  above,  the  calculated  response  of  the  Silver  Summit  Mine  drift  does  not 
agree  well  with  measurements  when  viscoelasticity  in  shear  is  assumed. 

As  to  the  finite  element  mesh  used  in  the  final  three-dimensional 
calculation,  some  auxiliary  calculations  have  shown  that  the  present  mesh  is 
definitely  too  coarse  to  calculate  accurately  the  extensometer  readings  and 
the  stress  concentrations  in  the  immediate  vicinity  of  the  cavity  in  the  early 
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stages  of  the  analysis.  The  mesh  is  adequate  for  determining  the  response 
(deflection)  of  the  back  and  the  floor  during  latter  stages  of  excavation.  It 
has  been  conclusively  shown  that  the  large  humps  present  in  the  calculated 
curves  around  the  18th  day  are  the  results  of  the  relative  coarseness  of  the 
mesh  over  the  back  in  the  early  stage  of  the  excavation.  When  Steps  4,  5, 
and  6  are  taken,  only  three  rows  of  elements  have  been  removed  so  that,  as 
seen  from  the  z-direction  (Figure  2-13),  there  are  only  three  elements  over 
the  back.  As  a  result,  the  mesh  over  the  back  is  too  coarse  during  this 
stage  of  the  calculation  to  represent  flexure  of  the  back  adequately. 

A  two-dimensional  analysis  along  the  longitudinal  axis  of  the  hoist 
room  Indicates  that  the  calculated  extensometer  readings  are  halved  by  using 
twice  as  many  elements  over  the  back.  The  present  writers  consider  that 
the  same  order  of  improvement  can  be  achieved  in  the  corresponding  three- 
dimensional  analysis  by  a  similar  refinement  in  mesh  over  the  back. 

To  illustrate  the  effects  of  including  the  changes  discussed  in  this 
section  on  the  results  of  analysis,  Figure  2-19  has  been  prepared.  The  shaded 
area  shows  approximately  the  extensometer  deflections  which  the  present  writers 
think  would  be  obtained.  The  reason  for  indicating  a  range  of  results  rather 
than  individual  curves  is  that  there  is  still  some  uncertainty  about  the 
elastic  stiffnesses  of  the  various  rock  types.  It  appears  that  the  values  of 
Young's  moduli,  listed  in  Table  2-2,  have  been  chosen  too  high. 
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RANGE  OF  COMPUTED  RESULTS  USING  IMPROVED  MODEL. 
EXACT  RESULT  DEPENDS  ON  ELASTIC  STIFFNESS 
ASSUMED  FOR  ROCK. 
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FIGURE  2-19.  AUTHORS'  ESTIMATE  OF  RESULTS  WHICH  WOULD  BE  OBTAINED  IF  MODEL 
WERE  TO  BE  REFINED  AND  A  SECOND  ANALYSIS  PERFORMED 
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2.5  SUMMARY  AND  CONCLUSIONS 

One  of  the  purposes  of  the  analysis,  which  involves  a  three- 
dimensional  viscoelastic  model  of  the  hoist  room  and  neighboring  rock,  is 
to  compare  the  measurements  with  calculations.  Another  purpose,  which  is 
considered  equally  important,  is  to  indicate  the  procedure  by  which  finite 
element  models  are  derived  from  geologic  maps,  from  laboratory  and  in  situ 
measurements  of  rock  properties,  and  from  the  judgment  of  geologists  and 
mining  engineers  acquainted  with  the  area.  Previous  comparisons  between 
measurements  and  analyses  which  have  been  made  by  Agbabian  Associates  and 
other  investigators  (Reference  2-8)  indicate  that  disagreement  may  exceed 
factors  of  3  or  4.  The  present  analysis  of  the  Caladay  Project  hoist  room 
should  be  evaluated  in  the  light  of  this  experience.  As  engineers  involved 
with  design  and  construction  of  mines  become  acquainted  with  the  requirements 
of  analysis,  information  will  be  developed  which  will  improve  the  quality  of 
future  analyses. 


There  are  several  ways  to  evaluate  the  present  study  of  the  Caladay 
Project  hoist  room.  One  way  is  to  compare  the  extensometer  measurements 
in  the  back  with  calculated  results,  as  is  done  in  Figure  2-15.  Another  is 
to  try  to  interpret  the  calculated  results  according  to  physical  intuition 
and  to  state  the  implications  of  this  type  of  analysis  for  design.  Such 
evaluation  is  limited  because  there  are  only  three  extensometer  measurements 
and  because  the  stiffnesses  of  the  rock  types  near  the  hoist  room  vary  so 
much  that  physical  intuition  is  not  very  helpful  in  explaining  the  outcome. 

The  analytic  model  also  has  some  shortcomings  which  make  comparison  difficult. 
These  are  considered  in  the  following  evaluation. 

The  present  analysis  of  the  Caladay  Project  hoist  room  involves 
the  following  steps : 

a.  Represent  the  hoist  room  geometry,  including  idealization  of 
continuously  variable  bedding  and  joint  planes,  into  three  or 
four  beds  with  average  or  representative  properties. 


48 


R-721 5-1 -2701 


b.  Determine  the  stress/strain  properties  of  each  rock  type  using 
laboratory  measurements,  in  situ  measurements,  and  judgment. 

The  argillite  is  assumed  to  be  represented  by  viscoelastic 
properties  because  this  is  the  simplest  mathematical  framework 
which  apparently  accounts  for  the  little  available  data. 

c.  Estimate  the  in  situ  stresses  based  on  judgment  and  a  few 
measurements  made  by  the  over-coring  technique.  In  situ 
stresses  are  idealized  as  being  uniform  on  the  outer  boundaries 
of  the  finite  element  mesh,  but  they  become  non-uniform  in  the 
interior  due  to  variation  in  the  stiffnesses  of  the  various  rock 
types  as  well  as  to  the  presence  of  the  hoist  room. 

Each  step  requires  some  idealizations  to  obtain  results  within  a  practical 
budget  and  some  assumptions  to  compensate  for  missing  data.  However,  these 
idealizations  and  assumptions  are  much  less  restrictive  than  those  which  have 
been  made  by  many  previous  investigators,  whose  work  is  often  limited  to  two- 
dimensional  or  three-dimensional  elastic  analyses. 

As  shown  in  Figure  2-15,  the  analysis  overestimates  two  of  the 
three  extensometer  measurements  by  a  factor  of  three  to  six,  depending  on  the 
stage  at  which  comparison  is  made.  The  third  extensometer  is  considered  to 
be  influenced  by  the  rope  raises  and  is  not  included  in  the  comparison. 

Factors  of  three  to  six  can  be  attributed  to  uncertainties  in  in  situ  moduli, 
Poisson's  ratio,  and  in  situ  stresses;  also,  the  measurements  contain  some 
uncertainty  which  has  not  been  evaluated.  The  fact  that  these  uncertainties 
can  be  overcome  to  the  extent  represented  by  Figure  2-15  is  moderately 
encourag i ng. 


Suggestions  as  to  how  to  improve  the  results  both  quantitatively 
and  qualitatively  are  made  in  Section  2.^*.  The  main  areas  in  which  the 
analysis  may  be  improved  are  application  of  in  situ  stresses,  definition  of 
viscoelastic  properties  and  refinement  of  the  finite  element  mesh.  These 
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suggestions  would  be  applied  to  subsequent  calculations  of  the  Caladay 
Project  hoist  room. 

a.  Replace  the  present  method  of  applying  boundary  conditions 
solely  by  external  forces  with  a  procedure  whereby  elements 
are  assigned  stresses.  The  assigned  stresses  would,  of  course, 
be  in  equilibrium  with  external  forces. 

b.  Reduce  the  zone  of  argillite  which  may  undergo  viscoelastic 
deformation  to  a  region  within  the  zone  if  influence  of  the 
hoist  room  excavation.  Consider  changing  the  viscoelastic 
model  such  that  there  is  time-dependent  extension  between  nodal 
points  corresponding  to  the  extensometers . 

c.  Refine  the  mesh  over  the  back  if  the  principal  objective  of  the 
calculation  is  to  obtain  accurate  deflections  of  the  back  at 
early  stages  of  the  analysis.  The  solution  at  later  stages 
appears  to  be  satisfactory. 

d.  Include  shotcrete  over  the  back,  which  acts  as  a  structural 
support  and  reduces  deflections. 

The  purpose  of  the  present  project  is  to  bring  the  most  advanced 
methods  of  structural  analysis  within  reach  of  mining  engineers  and  designers. 
The  analysis  of  the  Caladay  Project  hoist  room  indicates  how  far  the  present 
computer  program  will  go  toward  predicting  stresses,  roof  deflections,  and 
floor  heave  when  used  in  conjunction  with  the  kind  of  data  which  were  available 
for  this  study.  With  more  data,  the  analysis  would  be  better;  with  less  com¬ 
plete  data,  it  would  be  worse.  In  any  case,  there  will  be  some  important 
aspects  of  response,  such  as  stability,  which  must  be  assessed  by  noting  where 
stress,  strain  and  deflections  exceed  allowable  values.  The  designer  may  use 
the  computer  program  to  evaluate  several  candidate  designs  on  the  basis  of 
deflections  and  stress  concentration,  from  which  he  may  infer  that  rockfalls 
are  more  likely  with  one  design  than  another.  However,  the  final  decisions 
are  still  based  on  the  experience  of  the  designer,  of  which  calculations  such 
as  the  Caladay  Project  hoist  room  may  be  considered  a  part. 
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APPLICATION  OF  FINITE  ELEMENT  THEORY 

This  section  discusses  the  present  formulation  of  equations  of 
equilibrium.  The  provisions  to  extend  this  formulation  to  large  deformations  is 
also  described.  Then  the  types  of  elements,  including  truss,  beam  plane  strain 
and  axi symmet ri  c ,  three-di mens ional  and  thick  shell  elements  are  described.  In 
addition,  a  new  element  for  representing  slip  and  debonding  along  planar  joints 
is  described. 

3 . 1  SOLUTION  OF  NONLINEAR  EQUATIONS  OF  EQUILIBRIUM 

The  matrix  equation  of  equilibrium  for  a  structural  system  with 
material  nonlinearity  is: 


K  (u)  u  -  P  (3-D 

where  the  instantaneous  stiffness  matrix  (K)  is  a  nonlinear  function  of  the 
displacement  vector  (u) .  P  is  the  vector  containing  external  loads.  There 
exist  numerous  methods  of  solving  the  above  system  of  nonlinear  equations.  In 
general,  these  methods  can  be  divided  into  two  classes:  iterative  methods  and 
incremental  methods. 

Iterative  methods  apply  the  total  load  initially  and  approach  the 
solution  by  modifying  the  stiffness  matrix  and/or  modifying  the  load  vector. 
Modification  of  the  stiffness  matrix  in  general  accelerates  the  convergence 
but  is  computationally  costly.  The  optimum  may  be  achieved  by  occasional 
stiffness  reformulation. 
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In  incremental  methods  the  loads  are  appl ied  in  several  steps  and  an 
incremental  form  of  the  Equation  3_1  is  solved. 


K  uU  =  £P 

-n  --n+l  --n+1 


(3-2) 


where 


iu  ±l  =  U  ,  .  “  U 

--n+l  -n+l  -n 


£P  .  =  P  Al  -  P 

--n+l  -n+l  -n 


It  is  important  to  note  that  the  stiffness  matrix  K  can  only  be  formed  based 

n  ' 

on  the  displacement  vector  from  the  previous  step  u  which  creates  some 

**  n 

step-wise  error.  In  this  simple  incremental  technique  the  step-wise  errors 
can  accumulate  and  lead  to  considerable  total  error.  To  prevent  this  accumula¬ 
tion  of  the  step-wise  error,  a  modified  form  of  the  load  vector  is  used. 


K  Au 

-n  --n+l 


■  -n+i  -  -F„ 


(3-3) 


where 


=  Total  load  vector  at  the  end  of  the  (n+l)th  step 

=  Vector  of  the  internal  resisting  forces  at  the  end  of  the 
th 

n  s  tep 


By  using  this  method  of  load  vector  correction  the  equilibrium  is  satisfied  at 
the  beginning  of  each  incremental  step  and  thereby  the  accumulation  of  the 
step-wise  error  is  prevented.  Satisfaction  of  equilibrium  is  assured  in  spite 
of  errors  or  approximations  in  the  stiffness  matrix  and,  therefore,  the  reformu- 
lotion  of  the  stiffness  matrix  is  not  lequired  at  every  step.  However,  the 
error  in  each  step  is  directly  dependent  on  the  approximation  of  instantaneous 
stiffness  matrix. 
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An  alternative  method  is  to  apply  the  total  force  from  the  beginning, 
in  which  case  the  load  in  Equation  3-3  will  be 


P 

n 


total 


(n=l  , . . .  ,N) 


It  should  be  noted  that  the  application  of  the  total  loads  makes  this  method 
equivalent  to  an  iterative  scheme  with  load  vector  correction  which  was 
previously  discussed.  However,  the  loads  in  general  have  a  specified  history 
dictated  by  the  sequence  of  application,  sequence  of  construction  and  excava¬ 
tion,  and  the  time  phenomenon  associated  with  the  viscous  material  properties. 
In  most  practical  problems,  the  specified  history  of  loading  is  a  series  of 
step  functions.  This  is  true  in  case  of  construction  and  excavation  which  can 
be  considered  as  a  discontinuity  in  force-displacement  relation  and  an  abrupt 
change  in  the  instantaneous  stiffness  matrix. 

An  efficient  scheme  is  to  apply  the  total  of  the  step-wise  loading 
at  each  stage  and  then  carry  out  several  iterations  with  occassional  stiffness 
reformulation  to  accelerate  the  convergence.  This  scheme  is  summarized  in  the 
following  steps . 

For  each  s  tep : 

a.  Compute  u  =  u  .  +  Au  (for  first  step:  u  =  0) 

K  -n  -n-1  -n  o 

b.  Compute  the  strains  (e  )  or  strain  increments  (Ac  )  using  the 

-  n  -n 

derivatives  of  the  shape  functions  for  each  element  which  have 
been  initially  computed  and  stored 

c.  1.  For  time-independent  materials  compute  the  stress  (o  )  and 

the  i ns  tar taneous  s tress-s t ra i n  relations  (C  )  (see  section  on 

-n 

material  properties) 

2.  For  visco-elastic  elements 

£ 

i.  Compute  stresses  a  =  C  (c  -  e  ,)  where  C  is 
r  -n  -  -n  -n-l 

the  elastic  stress-strain  matrix,  c  is  the  total  strain 

c  n 

and  £n_j  is  the  total  creep  strain 
5*1 
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"•  Us!ng  the  Presses,  compute  ^  which  is  the  total 
creep  strain  at  the  end  of  the  time  step  (see  section  on 
material  properties) 

iii.  Compute  effective  stress 


o  =  C  (t  -  f  L) 

n  "  -n  -n 

d-  Compute  the  internal  resisting  forces  from  the  stresses  (effec¬ 
tive  stresses  for  viscoelastic  elements).  If  it  is  stiffness 
update  cycle,  compute  stiffness  matrix. 

e-  Solve  Equation  3-3  to  obtain  Aun+].  Compute  Y  =  ||  Aur+  || . 

f-  If  a  specified  number  of  iteration  has  been  reached  or  if 

Y  :  c  (e  is  a  specified  quantity),  go  to  Step  g;  otherwise, 
go  to  Step  a  and  repeat  the  iteration. 

9*  Apply  the  next  loading  step  and  go  to  Step  a. 

3,2  EQUATIONS  OF  EQUILIBRIUM  FOR  LARGE  DEFORMATIONS 

The  method  of  large  deformation  finite  element  analysis  to  be  used  in 
the  present  computer  program  was  initially  introduced  by  Shari f i  and  Yates, 
Reference  3_1 . 

The  matrix  equations  of  equilibrium  (or  motion)  are  derived  from  an 
incremental  virtual  work  expression  and  the  original  configuration  of  the  finite 
element  system  is  taken  as  the  reference  configuration.  This  choice  of  the 
reference  state  eliminates  the  need  for  updating  of  the  coordinates  of  the 
nodal  points  which  is  computationally  a  costly  operation. 
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The  incremental  virtual  work  expression  is 


f  (t  +  _t  )  6 All  dA  -  f  S..  5 Ac .  .dV  =  f 
Jh  n  n  n  o  ij  'J  °  Jy 


S.  .  6AH.  .dV 
ij  U  o 


+  /  AS . .  6Ae i j  dV 


(3_I0 


where 


u  =  Component  of  displacement  vector 


Au  =  Incremental  component  of  displacement  vector 
k 

t  =  Component  of  traction  vector 
k 

/\t  =  Incremental  component  of  traction  vector 

k 

S.j  =  Component  of  Piola  stress  tensor 

AS  =  Incremental  component  of  Piola  stress  tensor 

ij 

Ae  =  Linear  component  of  incremental  strain  tensor 

1  j 

AH  =  Nonlinear  component  of  incremental  strain  tensor 

'  j 

The  stresses  and  traction  are  referred  to  the  area  and  the  coordinates 
of  the  original  configuration: 


Ac.  .  =  Au,  .  +  Au  .  +  u,  .Au,  .  +  u,  .Auk  . 

1  j  1  1 J  J  . 1  1  K>J  K’J  K» 


(3-5a) 


AH.  .  =  Au,  .Au,  . 

ij  k, 1  k , j 


(3-5b) 


The  total  and  incremental  displacements  within  each  element  can  be 
expressed  in  terms  of  the  nodal  point  values  of  the  displacements  through 
shape  functions 
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where  H  is  the  vector  of  the  shape  functions  and  U.  and  AU  are  the 

-  1  — 1 

vectors  of  the  total  and  incremental  displacements  of  the  element  nodes. 

Without  loss  of  generality  the  remaining  part  of  this  section  will 
be  devoted  to  the  derivation  of  the  appropriate  matrices  for  the  two- 
dimensional  quadrilateral  element. 

Substituting  the  Equation  3-6  into  Equation  3~5  will  result  in 
the  following  expressions  for  the  strain  increments  in  terms  of  the  nodal 
point  displacements: 


-  ’y  Ayy  +  <fi'y  Yx>  (S'y  *  <B.y  V  (H.y  My) 

(3-7) 

B’x  4-y  +  -'y  A-x  +  (By  ’y  A-x^  +  <b-y  Wy) 

<H-X  A~y)  *  <S'x  ^  (S’y  AV  +  <%  A-y) 

The  above  equations  can  be  written  as  matrix  form 

Ac  =  B  (I  +  E)  AU  =  B  AU  (3-8) 

where  I  is  the  identity  matrix  and 
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§  is  the  usual  strain-displacement  matrix  for  infinitesimal  deforma¬ 
tion  and  E  is  the  large  deformation  contribution. 

The  linear  and  geometric  stiffness  matrices  imd  the  load  correction 

vector  are 

!5e  =  J  §T  C  B  dvm  (3-11) 

Vm 

bg  =  /  §T  §  §  dvm  (3-, 2) 

m 

f  =  /  §T  s  dvm  (3., 3) 

m 

where  v^  is  the  area  of  each  element  in  original  configuration,  and  C  is 
the  instantaneous  stress-strain  relation 

AS  =  C  A  c  (3_i4) 

Finally,  the  matrix  equation  of  equilibrium  is 

(Ke  +  Kg)  AU  =  R  -  F  (3-15) 
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It  is  important  to  note  that  for  the  computation  of  the  above  matrices 
only  the  derivatives  of  the  shape  functions  H,x  and  H,y  at  original  geometry 
of  each  element  is  required.  Therefore,  these  derivatives  at  integration 
points,  can  be  computed  in  the  first  part  of  the  program. 

3.3  STRUCTURAL  FINITE  ELEMENTS 

A  description  of  the  structural  elements  incorporated  in  this  computer 
program  are  given  here.  The  beam  and  thick  s he  1 1  elements  have  linear  elastic 
properties.  All  other  elements  are  capable  of  representing  nonlinear  properties 

3-3.1  THREE-DIMENSIONAL  TRUSS  ELEMENTS 

The  truss  element  is  the  conventional  space  truss  member  which  can 
resist  compression  or  tension  along  its  axis.  It  can  also  be  used  to  model 
bolts.  The  truss  member  is  subject  to  three  translations  at  each  end  of  the 
member  as  shown  in  Figure  3" 1  a .  The  member  stiffness  matrix  is  of  order 
6x6.  The  material  and  geometrical  properties  are  defined  by  the  tangent 
Young's  modulus,  and  the  cross-sectional  area  of  the  element. 

3.3.2  THREE-DIMENSIONAL  BEAM  ELEMENTS 

The  three-dimensional  beam  element  is  subject  to  three  translations 
and  three  rotations  at  each  end  of  the  member.  The  generalized  forces  and  the 
generalized  displacements  associated  with  the  s i x-degrees-of- freedom  (DOF)  at 
each  end  are  shown  in  Figure  3~ lb. 

The  geometrical  properties  of  the  beam  element  are  specified  by  an 
axial  and  two  shear  areas  and  three  principal  moments  of  inertia,  two  associated 
with  bending  and  one  with  torsion.  Young's  modulus  and  Poisson's  ratio  are 
required  to  define  the  material  properties  of  the  beam  element. 

The  element  stiffness  matrix  is  of  order  12  x  12  and  is  obtained 
from  the  classical  beam  theory  including  the  effects  of  the  shear  deformations. 
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AJA2I2I 


(a)  THE  THREE-DI  MENS  I  ONAL  TRUSS  ELEMENT 


AJA2122 


(b )  TIIC  THREE  -  D I  MENS  I  ONAL  13 LAM  E’  El'.EIIT. 
FIGURE  3-1.  TRUSS  AND  BEAM  ELEMENT 
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It  is  important  to  note  that  for  the  computation  of  the  above  matrices 
only  the  derivatives  of  the  shape  functions  H,x  and  H,y  at  original  geometry 
of  each  element  is  required.  Therefore,  these  derivatives  at  integration 
points,  can  be  computed  in  the  first  part  of  the  program. 

3.3  STRUCTURAL  FINITE  ELEMENTS 

A  description  of  the  structural  elements  incorporated  in  this  computer 
program  are  given  here.  The  beam  and  thick  shell  elements  have  linear  elastic 
properties.  All  other  elements  are  capable  of  representing  nonlinear  properties. 

3.3.1  THREE-DIMENSIONAL  TRUSS  ELEMENTS 

The  truss  element  is  the  conventional  space  truss  member  which  can 
res i s  t  comp  ression  or  tension  along  its  axis.  It  can  also  be  used  to  mode  1 
bolts.  The  truss  member  is  subject  to  three  translations  at  each  end  of  the 
member  as  shown  in  Figure  3~ 1  a .  The  member  stiffness  matrix  is  of  order 
6x6.  The  material  and  geometrical  properties  arc:  defined  by  the  tangent 
Young's  modulus,  and  the  cross-sectional  area  of  the  element. 

3.3.2  THREE-DIMENSIONAL  BEAM  ELEMENTS 

The  three-dimensional  beam  element  is  subject  to  three  translations 
and  three  rotations  at  each  end  of  the  member.  The  generalized  forces  and  the 
generalized  displacements  associated  with  the  s i x-degrees-o f- f reedom  (DOF)  at 
each  end  a?-e  shown  in  Figure  3-lb. 

The  geometrical  properties  of  the  beam  element  are  specified  by  an 
axial  and  two  shear  areas  and  three  principal  moments  of  inertia,  two  associated 
with  bending  and  one  with  torsion.  Young's  modulus  and  Poisson's  ratio  are 
required  to  define  the  material  properties  of  the  beam  element. 

The  element  stiffness  matrix  is  of  order  12  x  12  and  is  obtained 
from  the  classical  beam  theory  including  the  effects  of  the  shear  deformations. 
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AJA2I21 


(a)  THE  THREE-DIMENSIONAL  TRUSS  ELEMENT 


/ 


AJt.2122 

(b)  THE  THREE-DIMENSIONAL  REAM  ELEMENT. 


FIGURE  3-1.  TRUSS  AND  BEAM  ELEMENT 
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A  provision  for  the  member  end  boundary  conditions  accounts  for 
hinges  and  other  releases. 


3.3.3  TWO-DIMENSIONAL  PLANE  STRAIN  AND  AXI SYMMETRIC  ELEMENTS 

Quadrilateral  isoparametric  elements  will  be  used  in  the  computer 
program.  For  a  general  quadrilateral  element,  as  shown  in  Figure  3"2,  the 
local  and  global  coordinate  systems  are  related  by 


4 


£ 

h .  x. 

i  =1 

4 

£ 

i=l 

h.  y. 

(3-16) 


where  the  interpolation  functions  are  given  by 


^  =  1/4  (1-s)  (1-t) 
h2  =  1/4  (1+s)  (1-t) 
h3  =  1/4  (1+s)  (1+t) 
h4  =  1/4  (1-s)  (1+t) 

The  same  interpolation  functions  are  used  in  the  displacement 
approximation . 


ux  's,t^  =  1  N  uxi  +  hs  “1  +  h6  “2 
Uy  (s,t)  =  £  hj  uyj  ♦  h5  «3  +  h6 

where 

h5  =  (1-s2) 

h6  =  ('-t2) 


(3-18) 


hc  and  lv  are  the  incompatible  interpolation  functions. 
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For  two-dimensional  analysis  the  s t ra i n-d i sp 1 acemen t  equations  are 


3U* 

c  =  _ —  =  y  h  u  +  hr  a.  +  h-. 

XX  3x  1  ni,x  uxi  5,x  1  6,x  2 

3uv 

Eyy  =  “5y  '  E  h'  .y  uyi  +  h5,y  "3  +  hf>,v 


3u  3u 


f-  =  — *  +  — Y.  =  v  h,  u  .  +  Z  h.  u  .  +  Hr  a,  +  h/-  a-  +  hc  a_  + 

S<y  3y  3x  L  ni  ,y  uxi  i,x  uyi  5,y  1  6,y  2  5,x  3  6,; 


Or  Equation  3-19  can  be  written  in  matrix  form  as 


M,x 

0 

Ux 

e  =  B  U  = 

0 

H.y 

Uy 

H.y 

H.x 

—  — 

LT  "  J 

(3-19) 


(3-20) 


In  this  case  the  three  strains  are  related  to  the  eight  nodal  point  displace 
ments  and  four  coefficients  of  incompatible  displacement  functions  by  a 
3x12  matrix.  The  submatrices  in  Equation  3— 20  are  given  by 


-’x  =  Chl  .x  h2,x  h3,x  h4,x>h5,x’  h6,x] 

H.y  =  [h]>y  h2>y  h3>y  h4>y»  h5>y,  h6y] 


(3-21) 


The  element  stiffness  matrix  is  given  by  the  following  equation: 


k  =  /  BT  C  B  dv 

*v  - 


(3-22) 
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This  stiffness  matrix  which  is  12  x  12  is  reduced  to  8  x  8  by  elimination  of 
the  four  incompatible  modes  before  assembling  in  the  global  stiffness  matrix. 

3-3-4  THREE-DIMENSIONAL  ELEMENT 

For  an  arbitrary  eight-point  brick  element  shown  in  Figure  3“3,  the 
appropriate  displacement  approximations  are 

8 

\  ■  \i  *  h9“xl  +  h10°x2  *  hn  “x3 
8 

uy  ■  £  uyi  +  h9  “yl  f  h10  “y2  +  hu  “y3 
8 

Uz=  Uzi  +h9az1  +  h10az2  +  hn  az3 
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where 

h1  =  1/8  (1  +  a  (1  +  n)  (1  +  0 

h2  =  1/8  (1  -  K)  (1  +  n)  (1  +  0 

h3  =  1/8  (1  -  K)  (1  -  n)  0  +  0 

h4  =  1/8  (1  +  K)  (1  -  n)  0  +  0 

h5  =  1/8  (1  +  K)  (1  +  n)  (1  -  C)  (3_25) 

h6  =  1/8  (1  -  K)  (1  +  0)  (1  -  C) 

h?  =  1/8  (1  -  O  (1  -  n)  (1  -  C) 

h8  =  1/8  (1  +  O  (1  -  n)  (1  -  0 

hg  *  (1  -  C2) 

h10  =  0  “  o2) 

hn  =  (i  -  c2) 

The  first  eight  are  the  standard  compatible  interpolation  functions.  The  last 
three  are  incompatible  and  are  associated  with  linear  shear  and  normal  strains. 
The  nine  incompatible  modes  are  eliminated  at  the  element  stiffness  level  by 
static  condensation. 

3.3.5  THICK  SHELL  ELEMENTS 

The  thick  shell  element  described  here  was  initially  developed  by 
V/ilson,  et  a  1  .  ,  Reference  3~3- 

This  shell  element  is  a  16-node  curved  solid  element  shown  in 
Figure  3-4.  Each  node  has  three  unknown  displacements.  Therefore,  if  the 
shell  is  considered  as  a  two-dimensional  surface  there  are  six  unknowns  per 
point.  It  is  apparent  that  this  type  of  formulation  avoids  the  problems 
associated  with  the  sixth  degree  of  freedonr-the  normal  rotation  is  set  to 
zero  when  certain  finite  elements  are  used  in  the  idealization  of  shells. 
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IK 

here 

h,7  =  ?(i-  t2) 

h18  -  n  0  -  n2) 

h19  =  0  -  ?2)  (3-29) 

h20  =  En  0  -  S2) 
h21  *  nC  0  -  i2) 


The  motivation  for  addition  of  the  interpolation  functions  h,?  to 

is  to  increase  the  capability  of  the  element  in  producing  closer  approxi- 

sti ons  to  the  exact  displacements  under  simple  loadings,  thereby  increas.ng 

,e  convergence  to  exact  solution.  The  incompatible  Interpolation  functions 

-  .1  _  J _ rt^nrliiro  1  nrnmna t l b i 1 i t i GS  i H 

to 


;e  lu  eAULi  ^  ^ . - 

h  have  zero  values  at  the  nodes  and  produce  incompatibilities  in 
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isplacement  field  along  the  interelement  boundaries. 


,1|  JOINT  FINITE  ELEMENT 


The  joint  element  is  intended  to  represent  the  rock  joints,  faults, 
nterfaces  and  similar  discontinuities  in  continuum  systems.  The  joint 
dement  has  the  capability  of  representing  the  main  characteristics  of  the 
leforma t ion  behavior  of  the  rock  joints  such  as  debonding  and  slip.  The 
term  dehouding  means  the  ability  of  separation  of  the  two  blocks  of  continuum 
Kljacent  to  the  joint  surface  which  were  initially  in  contact.  Subsequent 
-ontact  can  also  develop  by  the  movement  of  the  two  blocks  towards  eac 
The  term  slip  means  the  relative  motion  along  the  joint  surface  or  fault  when 
the  shearing  force  exceeds  the  shear  strength  of  the  joint. 


Previous  attempts  have  been  made  to  develop  discrete  elements  to 
represent  the  joint  behavior.  Goodman,  Taylor  and  drekke  (Reference  3-2) 
developed  a  simple  rectangular,  two-dimensional  element  with  eight  degrees 
of  freedom.  This  element  has  no  thickness,  and  therefore  the  adjacent  blocks 
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of  continuum  elements  cun  penetrate  into  each  other.  Zienkiewicz,  et  al., 
(Reference  3“*0  advocate  the  use  of  continuum  isoparametric  elements  with  a 
simple  nonlinear  material  property  for  shear  and  normal  stresses,  assuming 
uniform  strain  in  the  thickness  direction.  Numerical  difficulties  may  arise 
from  ill  conditioning  of  the  stiffness  matrix  due  to  very  large  off-diagonal 


terms  or  very  small  diagonal  terms  which  are  generated  by  these  elements  in 
certain  cases. 

To  avoid  such  numerical  problems  a  new  joint  element  is  developed 
below,  which  uses  relative  displacements  as  the  independent  degrees  of  freedom. 
The  displacement  degrees  of  freedom  of  one  side  of  the  slip  surface  are 
transformed  into  the  relative  displacements  between  the  two  sides  of  the  slip 
surface.  The  transformation  relations  are  as  follows: 


The  superscripts  T  and  B  refer  to  the  top  and  bottom  elements 
with  respect  to  the  slip  surface  respectively.  As  shown  in  Figure  3"5,  those 
degrees  of  freedom  of  the  upper  element  which  are  on  the  slip  surface  are 
transformed  but  the  degrees  of  freedom  of  the  lower  element  are  the  original 
displacement  quantities. 
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The  substitution  of  Equations  3"30  and  3-31  into  Equation  3-32  results 
in  the  strain-displacement  relation  for  the  element 


(3-33) 


The  stresses  and  the  strains  are  related  through  the  following 
material  property  matrix  C. 


(3-34) 


g  -  C  e 

In  general  the  above  stress/strain  relationship  for  rock  joints  is 
nonlinear,  the  details  of  which  are  given  in  Section  k. 


The  stiffness  matrix  for  the  joint  element  is  formed  in  n-s 
coordinate  system; 

-ns  =  Jv  ^  §  dv  (3-35) 
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and  transformed  to  the  x-y  coordinate  system  as  follows 


k  =  T  k  T 
-  -ns  - 


(3-36) 


where  T  is  the  transformation  matrix  containing  the  direction  cosines. 


2 B 1 )  2(A3  +  B2)  (A ,  -  2B1 )  (A  +  B.,) 

2(A2  +  2B] )  (A  +  B  )  (A  +  2B  ) 


Symmet  r i c 


2 (A ,  -  2B] )  2(A3  +  B2) 

2(A2  +  2B2) 


(3-37) 


where 


2  2 

C  a  +  C  b 
ss  nn 

2  2 

C  b  +  C  a 
ss  nn 


(C  -  C  )  ab 
nn  ss 


■  r(,ij  -  xi> 

■  f  (yJ  '  yi> 


B.  =  C  ab 

1  ns 

B  =  C  (a2  -  b2) 

2  ns 
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SECTION  k 

REPRESENTATION  OF  PROPERTIES  OF  ROCK,  INCLUDING  ANISOTROPY, 
INELASTICITY,  RATE  EFFECTS  AND 
PROPERTIES  OF  FAULTS  OR  JOINTS 

The  first  part  of  this  section  describes  homogeneous  properties  of 
rock  which  are  available  in  the  AA  computer  program.  As  used  here,  "homogeneous 
refers  to  properties  which  can  reasonably  be  averaged  over  several  feet,  which 
is  a  typical  dimension  of  a  finite  element  in  applications  to  mining  engineer¬ 
ing.  Inhomogeneous  properties  of  rock  masses,  such  as  those  caused  by  faulting, 
are  treated  by  a  separate  procedure  which  is  described  in  the  second  part. 

The  topics  which  are  covered  below  include  the  following: 


a .  I ne I  as  t i c i ty 

1  .  Variable  modulus 

2.  Variable  modulus  with  perfect  plasticity 

3.  Variable  modulus  with  perfectly  plastic  fracture 
criterion  and  strain  hardening  cap 

b.  An i sotropy 

1.  Variable  modulus  with  anisotropic  fracture  criterion 

based  on  the  hypothesis  of  Jaeger  (plane  geometry  only) 

2.  Variable  modulus  with  anisotropic  yield  criterion  based 

on  the  hypothesis  of  Hill 

c.  Rate  Effects  (isotropic  Only) 

1.  Creep  (series  of  Kelvin  elements) 

2.  Viscoplasticity  (based  on  work  of  Perzyna) 
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d.  Joint  Properties 


1  .  Di 1  a  tan  t 


2 .  Nond i I i tan  t 


*4.1  HOMOGENEOUS  PROPERTIES 


These  models,  which  are  also  summarized  in  Table  *4-1,  afford  some 
material  descriptions  which  probably  cannot  be  effective!  used  at  present 
owing  to  the  lack  of  experimental  data  on  rock  from  particular  mines.  As 
analysis  increases  in  effectiveness,  it  is  expected  that  more  data  will  be 
routinely  obtained,  thus  enabling  more  sophisticated  material  models  to  be 
used.  Advantages  and  disadvantages  of  these  models  are  summarized  in 
Table  <4-2.  Typical  laboratory  properties  for  rocks  are  shown  in  Figure  *4-1 

TABLE  *4-1.  SUMMARY  OF  AVAILABLE  MATERIAL  PROPERTIES 


Type  of  Model 


Constant  elastic  moduli 

Variable  modul i 

Elastic/plastic  (fixed 
static  yield  surface) 

Constant  modul I 

Variable  modu! i 

Cap  model 


.tions  may  be  used  simultaneously,  but 
s  may  not  be  mixed,  except  where  noted. 
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TABLE  k-2.  ADVANTAGES  AND  DISADVANTAGES  OF  EACH  MODEL 


Advantages  Disadvantages 

A.  E I  as 1 1 c- I dea 1 )y  Plastic 


Simple  to  fit  to  data 
Approximates  most  features  of  data 
G  ■  Const.  G  =  G(P  )  and  associated 
flow  rule  theoretically  correct 


May  not  fit  all  available  data 
Cannot  match  ^riaxial  test 
Other  treatments  of  G  can  lead 
to  possible  paths  of  energy 
generati on 

For  non  assoc i ated  flow  rule  no 
general  uniqueness  theorem 


B .  Variable  Modu I i 


Best  fit  of  data 

Only  model  with  repeated  hysteresis 
within  failure  envelope 
Ideal  for  finite  element 
Computationally  simple 
Relatively  easy  to  f i t 

C. 


Satisfies  all  rigorous  theoretical 
requi remen ts 

Reasonably  good  fit  of  data 
Effective  control  of  dilatancy 


Restricted  to  near-proportional 
loading  (in  shear) 

For  nor,p ropor t  i onal  loading  paths 
no  uniqueness  theorem 
Additional  quantity  must  be  stored 
at  each  grid  point 

Cap  Model 

Indirect  approach  needed  to  fit 
data 

Relatively  complicated 
Additional  quantity,  the  strain 
hardening  parameter,  must  be 
stored  at  each  grid  point 


D.  Viscoelastic 


Simple  to  fit;  to  data  Requires  sophisticated  testing 

Approximates  features  to  define  viscous  coefficients 

of  data  for  some  rocks  for  multi-axial  loading. 

Does  not  account  for  deterioration  in 
s  trength  wi th  time 

E.  Viscoplastic 

Approximates  some  features  of  Requires  sophisticated  testing 

data  Including  shear  strength, 
stick-slip  phenomenon 
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(a) 

HYDROSTAT 
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failure  envelope 


TRIAXIAL  COMPRESSION 


TRIAXIAL  COMPRESSION 


AJA3751 


FIGURE  4-1.  TYPICAL  LABORATORY  OATA  ON  ROCK  FROM  WHICH  CONSTITUTIVE 
EQUATIONS  ARE  DERIVED  VE 
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The  programming  of  the  material  property  subroutines  is  arranged  to 
provide  maximum  flexibility  and  ease  in  changing  properties.  This  is  done  by 
performing  each  separate  function  in  the  materia!  properties  description  by  a 
separate  subroutine.  Thus  separate  subroutines  are  provided  for  the  following 
purposes : 


Computing  variable  moduli 

Computing  derivatives  of  the  yield  functions  with  respect 
to  i ts  argumen ts . 

Testing  for  yielding  or  fracture 

Adji.  ting  stresses  for  viscoelastic  relaxation. 

Adjusting  stresses  for  viscoplastic  relaxation. 

Transforming  strain  increments  to  principal  axes  of  orthotrophy 
and  transforming  the  matrix  of  stress/strain  coefficients  to 
the  global  axes. 

There  is  considerable  interdependence  between  the  inelastic  and  anisotropic 
capabilities.  Wherever  possible  this  interdependence  is  used  to  economize  on 
programming.  Logic  diagrams  for  each  subroutine  in  the  homogeneous  material 
property  package,  are  given  in  Appendix  A. 

*4.1.  I  INELASTICITY  FOR  ISOTROPIC  MATE  RIALS- -VAR  I  ABLE  MODULUS 

Inelasticity  in  isotropic  materials  is  represented  through  variable 

bulk  and  shear  moduli  and  through  plasticity  theory.  The  bulk  modulus  B 

is  assumed  to  depend  on  the  current  value  of  elastic  volumetric  strain  p  and 

its  previous  maximum  value  p 
r  max 

FOR  LOADING  (0  >  p  >  pm  ) 

Hlo  X 

B  =  (B  -  B  )  exp  (  — )  (4-1) 

m  o  \p,  ! 
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FUR  UNLOADING  OR  KLLOAD I NG  (0  <  p  <  Pmax) 

D  =  13  +  tB  -  B  )(  — ) 
u  in  u  \u  2 1 

whe  re 

(b  +  (b  -  b 

lo  m  o  \  U  / 

B  =  the  I esser  of  \ 

“ 

\  111 

FOR  LOADING  OR  UNLOADING/RELOADING  IN  TENSION  (u  ^  0) 


Application  of  this  model  to  a  granitic  rock  is  illustrated 
Figure  h-2.  Specific  parameters  for  this  rock  are 

,  6 

U  =  7 • 6  x  1 0  ps i 

m 

B  =  1 .205  x  106  psi 

o 

p,  -  0.0275 

u2  =  0.05 

The  shear  modulus  G  is  also  assumed  to  depend  on  p  and 


-1-2701 


(4-2) 


[b-r, 


P 

max 


FOR  LOADING  (0  <  p_  <  p) 


A 
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u  u 

(a)  MODEL  BULK  MODULUS  FOR  LOADING  (b)  MODEL  HYDROSTAT  FOR  LOADING 

AND  UNLOADING  AND  UNLOADING 


(c)  LOW  PRESSURE  BULK  MODULUS  COMPARED  (d)  MODEL  HUGONIOT  AND  HYDROSTAT 
WITH  DATA  (REFERENCE  4-1)  COMPARED  WITH  DATA 


FiGURE  4-2.  MODEL  BULK  MODULUS,  HYDROSTAT  AND  YIELD  POINT  IN  UNIAXIAL 
STRAIN  COMPARED  WITH  DATA 
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Theoretical  guidance  on  the  appropriate  functions  for  B  and  G 
is  provided  by  Walsh  (References  k-J ,  *t-8)  ,  who  postulates  that  the  effective 
modulus  differs  from  the  intrinsic  modulus  due  to  cracks  and  pores.  As  these 
are  closed  by  increasing  pressure,  the  eff  :ctive  modulus  tends  toward  the 
consolidated  value.  Walsh's  work  contains  parameters  which  are  not  retained 
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PRESENT  VARIABLE  G  HODEL 


SYMBOL  REFERENCE  TYPE  OF  ROCK 
A  \  UESTERLEY  GRANITE 

♦  {  ROCKPORT  GRANITE 

>  I4-S 

▼  I  STONE  MOUNTAIN  GRANITE 

■  )  PORTVILLE  GRANITE 

•  <t-6  NTS  GRANITE 


PRESSURE,  Kb 


SHEAR  MODULUS  VERSUS  PRESSURE  FOR  NTS  GRANITE 


— 1 - 

'1  I1  '  '-—l  - 

f 

■ 

▼ 

■ 

♦ 

V 

♦ 

■ 

♦ 

■ 

♦ 

A 

A 

♦ 

A 

PRESENT  CONSTANT  G  MOOEL 
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in  the  following  empirical  expressions  for  the  effective  bulk  modulus.  How¬ 
ever,  the  basic  concept  is  retained.  Also,  the  present  model  for  the  effective 
shear  modulus  merely  follows  Walsh's  concept.  The  idea  of  coupling  the  shear 
stiffness  to  volumetric  strain  is  proposed  in  Refererr.e  4-9  and  carries  the 
danger  that  energy  might  be  extracted  from  the  model  by  hydrostatic  com¬ 
press  icn,  followed  by  shearing,  followed  by  releasing  the  pressure  and  finally 
by  releasing  the  shear.  This  danger  is  avoided  by  assuming  that  friction 
prevents  cracks  from  reopening  during  unloading  So  that  the  largest  value  of 
G  reached  on  loading  is  retained  during  subsequent  uni oad i ng/rel oad i ng . 

Under  tnese  conditions  a  material  may  dissipate  energy  in  shear  duiing 
loading  and  unloading  cycles  but  can  never  produce  additional  energy. 

4.1.2  INELASTICITY  FOR  ISOTROPIC  MATERIALS— VARIABLE  MODULI  WITH  PLASTICITY 

The  present  adaptation  of  plasticity  theory  is  based  on  work  of 
References  4-9  through  4-13.  The  model  consists  of  a  yield  criterion 

f(a..,  L)  =  0  (4-8) 

where  Lisa  function  of  plastic  strain, 

and  a  plastic  flow  rule  in  which  f  is  regarded  as  a  potential  function 

dc?.  -  r<  (4-9) 

,j  30.. 

The  incremental  stress  is  related  to  the  elastic  component  of  incremental 
strain  by  Equation  4-7-  Defining 


/Il\ 

Substituting  Equations  *-9  and  *-10  into  Equation  *-7  leads  to 
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Mdekk ' A 


3f 


Do 


U 


U 


+  2G^de. 


3°ij/ 


('•-ID 


wlie  re 

X  =  li  -  j  G 

If  the  yield  criterion  is  satisfied,  the  stress  state  must  lie  on  the  surface 
defined  by  f  in  Equation  *-8.  The  mathematical  statement  of  this  constraint 
i  s 


df 


D  f  ,  ,  3  f  ,, 

t- —  do.  .  +  —  dL 
Do .  .  ij  DL 
'J 


0 


Substituting  Equation  *-11  into  *4-12  permits  solutions  for  A 


(*-12) 


A 


x(dckk)  fn+  2G  dciifii 
xfkkfu  +  2G  fijfij  +  R 


(*-13) 


where  R  is  a  strain-hardening  function  to  be  defined  below. 

Substitution  of  Equation  *-13  into  *-11  expresses  the  stress  increment  in 
terms  of  the  strain  increment. 


A 
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Specific  functional  forms  have  been  assumed  for  f.  These  contain 
empirical  constants  whose  values  can  be  selected  to  match  data  for  a  specific 
material.  The  forms  are 

Polynomial  1  in  0. . 


fl(°!J>  * 


V7!  •  £.  v' 

n=  1 

a5 


Pol ynomi a  I  2  i n  a . . 


f,(o.j)  = 


(V3!  -  {•.[' '  (’  *  4)]  *  j 

(V3!'  (ai  *  az> 


J1  >  b 


J,  <  b 


J 1  >  b 


J1  -  b 


(4-14) 


(4-15) 


Cap  (to  be  used  with  Polynomial  2) 


f2  '  (J,  -  V)2  +  p2(JJ  -  «)  '  0 


(4-16) 


in  which 


V  =  L  +  P  X (L)  X1  (L) 


=  [x(l)]2|i  +  PZ[X'(L)]2| 


(4-17) 


(4-18) 


x  (l)  = 


('-i)2] 


A  +  C 


L  <  B 


L  >  B 


(4-19a) 


89 


R-721 5- 1 -2701 


X  1  (L )  = 


j^C 

(° 


L  <  B 


L  >  B 


(4-i9b! 


The  hardening  parameter,  L,  is 


L  -  _f  g(j,,  V Jj)  dt 


(k-2G 


where 


xPI- 


(^-21 


/  \2  /  \2  /  v2 

fa  +  fa  +  fa 


(b-22 


The  hardening  parameter  R  is 


fa  -  Vs! 


3LU  \9a, 


(*»-23 


where 


Op  c^.  are  principal  stresses 


R  =  0  if  f  =  f . 


The  cap  parameters  and  stress/strain  relations  produced  by  the  cap 
model  are  shown  in  Figures  b-b  through  ^-6.  Data  on  strength  for  granite 
containing  various  degrees  of  cracking  are  shown  in  Figure  b~7 ,  which  illus¬ 
trate  the  adequacy  of  the  assumed  fracture  criteria. 


M 


STRESS,  KSI 
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1% 


RADIAL  STRESS,  KSI 


FIGURE  A-5.  HYDROSTATIC  AND  UNIAXIAL  STRAIN  BEHAVIOR  CAP  MODEL  FIT  FOR 
r.DAuiTir  matfria; 


STRESS  DIFFERENCE,  KSI 
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Ttie  incremental  stress/strain  equations  are  expressed  in  matrix  form 
as  foil ows : 

{da}  =  [ C ]  { d c }  (*>-24) 

where  dt_  is  the  total  increment  of  strain.  The  C  matrix  thus  contains 
generalized  tangent  moduli  and  can  be  used  in  forming  the  element  stiffness 
matrices.  For  the  models  described  above,  the  C  matrix  is  as  follows: 


95 
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The  C  matrix  is  clearly  composed  of  elastic  and  inelastic  parts 

C  =  ce  -  CP  (4-26) 

The  C  and  matrices  are  computed  separately.  The  reason  for  this  is 

efficiency  in  treating  isotropic  and  anisotropic  materials  with  the  same 
Fortran  statements. 

In  weakly  nonlinear  problems  it  is  possible  to  avoid  time-consuming 
reformulation  of  the  stiffness  matrix  by  introducing  nonlinearity  through 
the  load  vector.  The  method  is  an  extension  of  Lhe  equilibrium  equations 
given  in  Reference  b-29.  In  the  following  equation,  time  is  used  as  a 
parameter.  Number  of  load  step  could  be  used  instead. 

At  time  t,  the  total  change  in  complimentary  strain  energy  is 
equal  to  the  change  in  complimentary  work  done  by  nodal  point  forces. 

T=t  r  i=c 

X)  I  <e>T  (da}  dV  =  Yj  <u>t  ( dP }  (4-27) 

T=0  *vO I  T  T  t»0  T 


where 

<c>, 

{da} 

=  Element  strain  and  stress  increment 

<u>, 

{dP} 

=  Nodal  displacement  and  force  vector 

V 

=  Volume  of  finite  element 

T 

=  Arbitrary  instant  of  time 

The  strain/displacement  relation  is 


A 
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The  stress  increment  in  an  elastic/plastic  material  may  be  expressed  b> 
rewriting  Equation  as  follows. 


{do}  =  [C ] { de }  -  [da  } 

P 


(*»-29) 


elastic  correction 


where  now  C  =  C 


Substituting  Equations  2-28  and  2-29  into  Equation  2-27, 


5^  f  <u>  [B]T([C] [B] {du}  -  {da  }  IdV  =  £  <u>  {dP}  (A~30) 

t=0  •'vol  T  *  PT  t=0  T  T 


Noting  that  <U>T  may  eliminated  from  both  sides  and  that 


[ B ]  [C]  [B]  dV  =  [K] 


(4-30 


where  [ K ]  is  the  elastic  stiffness  matrix,  Equation  2-30  may  be  rewritten  as 


i=t-At 

£ 

1=0 


^[K]{dulT  -  J  [B]  T{da_i)^  dV^ 


C  T  T=t 

[l<]{du}Ai  -  I  [B]  | da  |  dV  =  Y.  {p}T  (l*-32) 


where  { d u }  A  ,  jda  |  equal  change  in  u  and  a  during  the  interval 
At  P  At  P 

t  -  A  t  to  t . 

Performing  the  indicated  summation,  assuming  stress  to  be  constant 
throughout  the  elenier*-  and  defining 


E  =  Y,  [ K] { du } 

i=:0 


(4-33) 
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A 

we  have 


lKHdu)At  -  1b!T|%I  V  -  (P)t  -  (E)  *  [B]T{0  \  V  tt-3'1) 

At  1  P;t-At 

For  small  increments  of  stress  and  time  in  a  step-by-step  integration,  the 
second  term  on  the  left  hand  side  is  neglected. 

The  expression  which  is  used  for  {E}  in  the  computer  program  is 
derived  as  follows.  The  recoverable  work  done  on  an  element  is  equated  to 
the  elastic  strain  energy  stored  in  the  element  by  the  following  equation. 


I<u>{E}  =  y  <e>[c]{e}V 


(A-35) 


or 


(E)  =  [BT] [C] {e}V 


(*•-36) 


Thus  Equation  4-34  may  be  rewritten 


[K] {du} 


At 


=  -  [BT]  j[CHe) 


(*♦-37) 


Following  Equation  4-29 


do . 


>■( 


Xdekk6ij  +  2G  dE 


"  (A  ~~  6.  .  +  2G 

IJ/  *  3okk  ij  9a../ 


A  (4-38) 


or 


da.  . 

U 


Ade, . 6. .  +  2G  de. .  -  doP. 

kk  i  j  i  j  i  j 


(4-39) 
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where 


dJ.  .  = 

|J 


3a 


=  0 


-f  6 .  .  +  2G 

ij  3a.. 


kk 


Ade 


3f 


kk  3a 


+  2G  de 


3f 


JU 


ij  3a. 


3f 


3  f 


3a 


kk 


3a 


+  2G 


3f 


LL 


3f 


n 


3a . 


•J 


3a .  . 
'J 


if  f  =  0 


(4-40) 


if  f  <  0 


After  each  integration  step,  Au.  is  given;  dt..  is  found. 

Based  on  de..  and  stress  at  previous  time,  f  is  checked 
and  da.j  is  calculated  according  to  whether  material  is 
elastic  or  plastic  by  adding  contribution  for  elastic  part 
to  contribution  from  plastic  part. 

Stiffness  is  always  based  on  elastic,  parameters,  i.e.,  A,  G. 
Plasticity  is  introduced  through  updating  of  stress  increment. 
Hence,  there  is  need  to  update  stiffness  matrix. 


INELASTICITY  FOR  ANISOTROPIC  MATERIALS— THEORY  OF  JAEGER 

The  fracture  of  anisotropic  rocks  is  the  subject  of  several  failure 
theories.  The  Walsh-Brace  theory  (Reference  4-l4)  assumes  that  failure  is 
tensile  in  nature  and  that  it  is  influenced  by  the  presence  of  preexisting 
cracks.  Some  of  the  cracks  are  assumed  to  be  small  and  randomly  oriented 
while  others  are  long  and  have  preferred  directions.  Extension  of  the  cracks 
is  postulated  to  occur  when  the  Griffith  criterion  (Reference  ^-15),  as 
modified  by  McClintock  and  Walsh  (Reference  *.-16)  to  account  for  friction  on 
the  crack  faces, is  satisfied. 


In  contras';,  Jaeger  (Reference  4-17)  assumes  the  material  to  fail  in 
shear  either  along  a  single  plane  of  weakness  or  within  the  matrix  material 
according  to  a  Mohr-Coulomb  type  of  failure  criterion  of  the  form 


T 


=  a 

o 


-  aa , 


(4-38) 


[il»j 


#ll\ 
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where 

t 

a 


Sheer  stress  on  plane  of  fracture 
Normal  stress  on  plane  of  fracture 
Cohesion,  angle  of  friction 


Th  heory  is  expressed  in  terms  of  principal  stresses  by  an  axis  rotation 
as  : o I  lows 


2ao  ’  2a1CT3 

°1  "  °3  =  - —  -■  (4-39) 

a,  -  yifTT 

Improved  agreement  with  experiment  is  obtained  i f  a^  is  assumed  to  vary 

ao  =  a2  '  a3  CCS  (2^  '  2')  (4-4°) 


where 


C  =  Counterclockwise  angle  from  the  direction  of  the  major  prin¬ 
cipal  stress  (<jj)  to  the  direction  of  the  bedding  planes 

=  The  orientation  of  3  at  which  a  is  minimum.  Usually 
assumed  to  be  equal  to  30  deg  ° 

The  angle  3  is  shown  in  Figure  4-8  along  with  angles  relating  it  to  global 
directions  in  the  finite  element  formulation.  McLamore  and  Grey  (Refer¬ 
ence  4-18)  have  obtained  satisfactory  agreement  between  experimental 
data  on  strengths  of  slates  and  shales  using  a  modification  of  Equation  4-25 
as  fol lows : 


Some  of  their  resul ts 
Figure  4-9. 


-  a^  [cos  2(5  "  3)  ] n 

-  a,,  [co  2(?  -  3)  ]n 

and  those  of  Brace  and 


for  0  <  5  <  3 

for  3  <  5  <  90° 

Walsh  are  illustrated  in 


(4-41) 


0 
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x  "  v  -  z  =  GLOBAL  AXES 
1  -  2  -  3  =  PRINCIPAL  STRESS  AXES 


FIGURE  k-8.  ORIENTATION  OF  PLANES  OF  WEAKNESS  DEFINING  ANISOTROPIC  BEHAVIOR 
ACCORDING  TO  JAEGER 
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COMPARISON  OF  VARIOUS  FAILURE  THEORIES 
WITH  SLATE  DATA  FOR  CONFINING  PRESSURES 
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FIGURE  4-9.  COMPARISON  BETWEEN  FAILURE  THEORIES  AND  EXPERIMENTAL  DATA 
FOR  ANISOTROPIC  ROCK  {REFERENCE  2-18) 
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The  present  application  of  the  hypotheses  in  Reference  4-18  is 
limited  to  plane  geometry.  It  does  not  take  into  account  the  effect  of  the 
intermediate  principal  stress  which  is  shown  in  Reference  4-19  to  play  an 
important  role  in  fracture  of  some  types  of  rocks.  More  work  is  needed  to 
remove  these  restrictions. 

The  first  step  is  to  determine  the  magnitude  of  the  principal 
stresses  ^ and  their  direction  as  specified  by  0,  the  counter¬ 
clockwise  angle  from  the  +X  global'  axis. 


o 


1,2 


o  + 
XX 


2 


± 


(4-42) 


0 


=  i-  Arctan  (-——2X— 
2  la  -  o 

\  y  * 


(4-43) 


Bilinear  stress/strain  relations  are  specified  by  the  user  in  terms 
of  Young's  moduli  and  Poisson's  ratio  in  directions  parallel  and  perpendicular 
to  the  bedding  planes  (oriented  at  a  relative  to  global  axes).  Thus, 
experimental  data  is  required  from  specimens  cut  orthogonal  to  bedding  plane 
and  at  angles  other  than  90°.  The  computer  program  transforms  the  various 
E  and  v  to  the  principal  directions  of  stress  and  modifies  them  to  account 
for  fracture.  These  parameters,  E]t  v]2>  v^,  E2,  \>2},  \>2  ,  etc.,  are 
assembled  into  a  matrix  relating  incremental  stress  and  strain  in  principal 
stress  axes.  The  relationship  between  incremental  stress  and  incremental 
strain  expressed  in  the  principal  axes  of  anisotropy  (principal  stress  axes) 
is  shown  in  Equation  2-24  where  C  is  given  by  Equation  4-44.  The  matrix  is 
then  transformed  through  the  angle  0  into  global  coordinates  for  inclusion 
in  the  element  stiffness  matrix.  An  illustration  of  the  bilinear  Young's 
modulus  approach  is  superposed  on  data  in  Figure  4-10. 
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IVES  FOR  GREEN  RIVER  SHALE-2  FOR  VARIOUS 
IES,  B  =  10  DEG  (REFERENCE  2-18) 
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where 


1  "  v23 V3 2 


V1 2V2 1 


'  V13V31 


v12v23v3i 


v13v21v32 


4.1.4  INELASTICITY  FOR  ANISOTROPIC  ROCK— THEORY  OF  HILL 

Hill  (Reference  4-9)  has  proposed  a  yield  criterion  for  anisotropic 
materials  whose  form  is  compatible  with  the  isotropic  plasticity  theory 
described  above.  If  the  stress  components  are  expressed  in  the  principal 
axes  of  anisotropy  (not  necessarily  principal  stress  axes),  the  yield  criterion 

can  be  expressed  by 
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in  which  it  is  assumed  that  the  A,  n,  4  axis  system  coincides  with  the  principal 
direction  of  anisotropy.  f  may  be  used  as  the  potential  function  described 
above.  Adaptation  of  this  theory  to  rock  is  described  in  Reference  4-20. 

The  elastic  behavior  of  the  material  may  be  prescribed  to  be  either 
isotropic  or  anisotropic.  If  it  is  isotropic,  the  quantities  B  and  G  may 
be  used.  If  it  is  anisotropic,  Young's  moduli  and  Poisson  ratios  E^,  v^. 
v^,  etc.  are  specified  in  the  principal  directions  of  anisotropy.  The  C 
matrix  (Equation  4-24)  which  relates  incremental  stress  to  incremental  strain 
is  thus  initially  expressed  in  the  principal  axes  of  anisotropy  and  is  sub¬ 
sequently  transformed  to  global  directions  of  the  finite  element  mesh. 

A.  1.5  RATE  EFFECTS— VISCOPLASTICITY 

This  method  of  incorporating  rate  sensitivity  equations  is  based  on 
Perzyna's  elastic-viscoplastic  model  (Reference  4-21)  which  is  a  generaliza¬ 
tion  of  an  earlier  model  proposed  by  Hohenemser  and  Prager  (Reference  4-22) . 

An  adaptation  of  the  cap  model  described  above  for  viscoplasticity  is  described 
in  Reference  4-23.  The  present  model  is  taken  from  Reference  4-24. 

A  linear  elastic,  rate  independent  region  is  bounded  by  a  static 
yield  cr i ter  ion 


f(J,,  JJ)  <  0  (4-46) 

within  which  Hooke's  Law  applies.  If  the  static  yield  criterion  is  satisfied 
or  exceeded 

f  (J, ,  Jp  >  0  (4-47) 

A  viscoplastic  strain  rate  is  assumed  to  develop  according  to  the  following 
flow  rule. 

Jij  ■  ^(f)  dr-  i*-^) 

'  J 
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where 


♦  (f) 

F 
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A  function  of  the  static  yield  criterion  f 


Assumed  potential 
flow  rule  is  used 
i s  performed  wi th 


function.  Presently,  a  nonassociated 
in  which  F  =  and  differentiation 

respect  to  the  stress  deviators  a!. 

'  J 


Empirical  viscoplastic  parameter 


It  is  assumed  in  the  present  work  that 


♦  (f) 


(-A- 


1 . 


1 . 


for 


f 


,n-1 


fJ2  -  S  '  2  0  J,  >  b 


n=1 


y[^2 


J2  •  a5  -  0 


J,  <  b 


(4-49) 


(4-50) 


If  f  <  0,  elastic  inviscid  stress/strain  relations  apply. 


Making  use  of  Hooke's  law 


a  =  Xe,  .  5 .  .  +  2G  c? . 

kk  i j  i  j 


(4-50 


and  expressing  the  elastic  deviatoric  strain  rate  by 
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/ll\ 

for  the  case  where  f  >  0,  the  stress  raie  may  be  expressed  by  the  followli 
equa  t i on  . 


Ac  6. .  +  2G|t . .  -  y 


T2-* —  -  il-y-  j  >  b 

25  vf"  )VJi 


Ac. .  6.  .  +  2G  E..  ^  -  1 

kk  ij  tj  I  a 


5  & 


(4-53) 


J 1  <  b 


If  time  is  considered  the  integration  parameter  and  the  time  step  is  At, 
the  incremental  stress  is 


Ade.,6..  +  2G  de..  -  y - 

kk  ij  ij  'U  +  a 1 J 1 


(4-54) 


The  absence  of  plastic  volume  strain  is  due  to  choosing  a  nonassoc ia ted  f  1 1 
rule  in  which  J.  does  not  appear  in  the  plastic  potential. 


Some  example  calculations,  which  are  summarized  in  Tables  4-3  and  4-4, 
are  shown  in  Figures  4-11  through  4-13*  Comparison  with  some  data  for  a 
tonal i te  is  shown  in  Figure  4- 1 4 .  The  comparison  illustrates  the  ability  of 
the  model  to  represent  increase  in  strength  with  strain  rate  and  its  inability 
to  represent  nonlinear  behavior  prior  to  yielding.  It  can  be  shown  that  the 
effects  of  viscoplasticity  can  be  accounted  for  entirely  by  a  correction  to  the 
load  vector  (Reference  4-24).  The  technique  is  similar  to  that  described  above 
for  rate- independent  plasticity  and  to  that  described  below  for  viscoelasticity. 
An  important  consequence  of  this  is  that  time-consuming  reformulation  of  the 
stiffness  matrix  is  avoided. 


[•] 


cr\  1-A 
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TABLE  A-3.  SUMMARY  OF  EXAMPLE  CALCULATIONS 


Case 

Type  of  Loading 

Stress  Ratio, 
°3/ol 

Strain  Rate,* 
i n. / i n . -sec 

1 

Proportional  loading 

0.293 

200.00 

2 

Proportional  loading 

0.293 

2.00 

3 

Proportional  loading 

0.293 

0.20 

A 

Proportional  loading 

0.293 

0.02 

5 

Proportional  loading 

0.200 

0.20 

6 

Proportional  loading 

0.  100 

0.20 

7 

Uniaxial  strain 

0.293 

0.02 

8 

Uniaxial  strain 

0.293 

0.20 

9 

Uniaxial  strain 

0.293 

2.00 

10 

Uniaxial  strain 

0.293 

200.00 

■ : F o i  proportional  loading,  these  are  elastic  strain  rates. 


TABLE  h-k.  PROPERTIES  USED  IN  PRESENT  EXAMPLES 


P  ropert i es 


Bu 1 k  modu 1  us 

B 

1.5  x  10 

Shear  modulus 

G 

10^  psi 

Cohes i on 

a  = 

o 

1 A50 

Angle  of  friction 

a,  = 

-0.1 

Viscoplastic  coefficient 

Y 

1  .0 

no 
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1.6  RATE  EFFECTS— VISCOELASTICITY 

The  total  strain  is  defined  to  be  the  sum  of  instantaneous  elastic 
and  viscoelastic  parts.  The  strain  is  further  divided  into  shear  and  volumetric 
components,  which  are  treated  separately  as  follows: 


f ij  f i j  "  ^ i j  Ekk 


kk 


~  (£kk)e  +  KJ 


kk 


(^-55a) 

(4-55b) 


In  the  computer  program  the  user  may  choose  to  have  either  elastic  or  visco¬ 
elastic  shear  deformation  and  to  have  either  elastic  or  viscoelastic  volumetric 


de format  ion . 


Kelvin,  Maxwell  and  three-parameter  fluid  models  are  available  as 
shown  in  Figure  1 5 .  To  simplify  the  following  discussion,  no  distinction  is 
made  between  volumetric  and  shear  components.  The  creep  rate  and  creep  strain 
at  t  +  At  may  be  expressed  as  follows: 


KELVIN 


•c  ,  c 
E  +  a1  E  =  a2o 


(^4- 56a  ) 


t  +  At  Et  exP  ( -a  1  A  t )  +  ot  —(I  -  exp  (-a  1  At) ) 


(^ - 56b ) 


where  ^  and  a,,  =  -  ,  E  and  n  are  the  spring  and  dashpot  constants 

respectively. 


MAXWELL 


c  =  a jO  +  a2a 


t  +  At 


(a0 At  +  a, )  a 
2  1  t  +  At 


a1°t 


(^-57a) 

(z+"57b) 


where  a.  =  rL  and  a.  =  — 
1  t  2  n 
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KELVIN 


-VW - MAXWELL 


3  PARAMETER 
FLUID 


15.  VISCOUS  MODELS  AVAILABLE  IN  THE 
PRESENT  COMPUTER  PROGRAM 
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THREE-PARAMETER  FLUID 


The  creep  strain 

.  1c  ,2c 
components  c  and  e 

and  dashpot 


of  the  three-parameter  fluid  is  divided  into  two 
which  represent  the  creep  strain  in  a  Kelvin  model 


c 


E 


1  C 
e  + 


2  c 
e 


(A-58a) 


1  c 

E  t  +  At 


1c  a2 

Et  exp (-a j At)  +  ot  —  (l  -  exp(-a1  At)) 


(A-58b) 


2c  2c  At 

f  t  +  At  f  t  +  At  a20t 


(4-58c) 


These  models  are  illustrated  in  Figure  4-15.  Various  aspects  of  their  proper¬ 
ties  are  discussed  in  Reference  4-35.  Equations  4~56a,  -57a,  and  -58a  are  in 
suitable  form  for  application  to  a  time-marching  integration  procedure. 


Typical  creep  and  creep  recovery  data  for  several  rocks  are  shown  in 
Figure  4-16.  Application  of  a  Kelvin  model  to  one  of  the  rocks  is  shown  in 
Figure  4-17.  By  suitably  varying  the  parameters  a1  and  a2>  a  range  of 
behavior  can  be  reproduced. 


The  finite  element  adaptation  of  these  models  uses  the  initial 
strain  approach  to  writing  the  equilibrium  equations  based  on  the  change  in 
internal  energy. 


U 


oTe 


dV 


J*  PTu  dS 


(4-59) 


The  s t re s s/ s t ra i n  relations  are  in  terms  of  a  matrix  C  of  elastic  stress/ 
strain  coefficients 


(aj 


[C] !c  +  dc  -  cH 


(4-60) 


Strain,  <  x  10' 


m\ 
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I  t  (diipUceiwnl) 

J.2r  («•) 
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(a)  CREEP  AND  CREEP  RECOVERY  (SLATE) —REFERENCE  4-28 


10  15 

Time,  sec  x  103 


6  hours 


20  25 


(b)  CREEP  UNDER  CONSTANT  STRESS  AND  RECOVERY  CURVES— REFERENCE  4-}0 
(AFTER  REFERENCE  4-31) 


FIGURE  4-16.  EXPERIMENTAL  DATA  ON  VISCOELASTIC  PROPERTIES  OF  SEVERAL  ROCKS 
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,  E,  =  ^76.2 
a)  1 

)  n,  =  3-950 


{  ' 
(n  =  1/2  i 


TIME,  10J  SEC 


FIGURE  4-17.  COMPARISON  BETWEEN  PRESENT  CREEP  MODEL  AND  EXPERIMENTAL  DATA 
(REFERENCES  4-30,  4-31) 
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where 


=  Stress  at  t 


Total  strain  (including  creep)  at  t  =  t 
=  Strain  increment  (including  creep)  during  At  =  t  -  t 
=  Accumulated  creep  to  t  =  t 


An  error  is  committed  in  Equation  A- 60  in  that  de  contains  both  elastic  and 
viscous  components,  whereas  it  is  treated  as  if  it  were  entirely  elastic. 
Defining  the  strain/displacement  transformation  matrix  B  as 


{ e }  =  [  B  ]  { u } 


(*-61) 


The  strain  energy  for  an  element  is 


U  =  duT  I  BTCBdV 


U  +  J  .  (Eo  -  cc) 


CBdV  u 


(*-62) 


Defining 


B  CBdV  =  k 


(*-63) 


where  K  is  the  element  stiffness  matrix  and  performing  a  variation  with 
respect  to  the  generalized  displacements  u  results  in 


kdu  =  P  -  F 


(*-6*) 


where 


(e  -  e  )T  CBdV 
6  U  c/ 
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It  IS  significant  that  the  element  sti  ffness*  matrix  consists  entirely  of 
linear  elastic  terms.  Thus,  the  assemblage  stiffness  matrix  needs  to  be 
formulated  only  once,  resulting  in  great  economy  of  computing.  It  is  neces¬ 
sary  to  store  the  components  of  the  creep  strain  at  t  for  use  in  the  next 
step.  ° 


This  approach  has  been  adapted  to  finite  element  for  rock  and 
concrete  (Reference  4-26).  Some  experimental  work  which  has  been  performed 
on  rock  is  compared  in  Reference  4-27  with  a  series  of  Kelvin  models. 

Reference  4-28  discusses  methods  of  accounting  separately  for  volumetric  and 
deviatoric  creep  strains. 

4.2  MATERIAL  PROPERTIES  OF  JOINTS 

This  section  describes  the  properties  which  may  be  assigned  to  joints. 
These  properties  consist  of  the  shearing  and  normal  stiffnesses  of  the  joints. 
They  correspond  physically  to  the  stiffness  and  strength  of  fault  gouge,  to 
the  roughness  of  the  joints  and  to  the  angles  of  slip  surfaces  relative  to  the 
principal  plane  of  the  joint.  They  are  classed  as  dilatant  if  shearing  pro¬ 
duces  joint  expansion  or  contraction  or  nondilatant  if  shearing  and  normal 
displacement  are  uncoupled.  The  properties  are  specified  by  the  user  in 
natural  coordinates  which  may  be  directions  parallel  and  perpendicular  either 
to  slip  surfaces  or  to  the  principal  plane  of  the  joint.  In  either  case, 
transformation  to  global  directions  is  automatically  performed  by  thr  program. 

As  in  the  case  of  homogeneous  material  properties,  the  joint  prop¬ 
erties  are  controlled  by  subroutines  CONECT  and  ELPL  through  modular  subroutines. 
Presently,  these  contain  built-in  joint  properties.  As  more  data  become 
available  on  properties  of  joints  and  the  present  joint  material  properties 
become  obsolete,  the  present  model  may  easily  be  modified  without  disturbing 
the  main  program  or  any  of  the  other  material  properties. 
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4.2.1  NONDILATANT  JOINTS 

This  class  of  joints  is  the  simplest  to  model  mathemat i :al ly  since 
there  is  no  volume  change  due  to  shearing  strains,  and  therefore  the  shear 
and  the  normal  components  of  deformation  are  uncoupled  and  the  stress-strain 
relations  are  as  follows: 


However,  C  and  C  are  nonlinear  functions, 
nn  ss 

In  stress-deformation  relationship  in  normal  direction,  three 
distinct  stages  can  be  recognized  (see  Figure  4-18); 

a.  Separation,  C  =  C  =0  when  c  i  0 

r  nn  ss  n 

b.  Crushing  of  the  surface  irregularities  or  the  compression  of 

the  material  in  the  fault  or  joint,  if  any  C  =  E 

J  nn  c 

(eC  <  e  <  0) .  For  smooth  surfaces  this  case  does  not  exist, 
n  n  c 

therefore  e  =0. 

n 

c.  Contact,  C  =  Ec.  (e  <  eC) 

nn  f  n  n 

It  is  important  to  note  that  very  high  values  can  be  assigned  to  E^  without 
any  numerical  problems  with  the  special  joint  element  described  in  Section  3- 

The  tangential  s t res s-s t ra i n  relationship  is  assumed  to  be  elastic- 
perfectly  plastic  using  a  Mohr-Coulomb  yield  criterion: 

C  =  G  o<c  +  o  tan  4? 

ss  s  n 

C  =  0  o  =  c .  +  o  tan  4> 

ss  s  n 

where  c  and  4>  are  the  cohesion  or  the  angle  of  friction. 


(4-66) 
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4.2.2  DILATANT  JOINTS 

Dilatancy  of  rock  joints  and  faults  are  very  complex  to  model 
mathematically;  however,  to  include  a  measure  of  dilatancy,  the  procedure 
developed  by  Goodman,  Dubois,  and  Brekke  (Reference  4-36)  is  used  here. 

Further  data  are  available  from  Reference  4-37. 

It  is  assumed  that  the  deformation  in  p-q  coordinate  system  shown 
in  Figure  4-19  is  nondilatant.  The  angle  y  between  the  p  direction  and  the 
joint  surface  is,  therefore,  defined  to  be  a  material  property  of  the  joint. 

The  stress-strain  relation  in  n-s  coordinate  is: 


(C  C2  +  C  s2) 
qq  pp 


(c  -  c  )  sc 
qq  pp 


(c  -  c  )  sc 
qq  pp 


(c  s2  +  c  c2)| 
qq  pp 


(4-67) 


where 

C  =  cos  y 
S  -  Sin  y 


4.2.3  DEMONSTRATION  OF  JOINT  ELEMENT 

The  example  problem  in  Figure  4-20  demonstrates  the  behavior  of  the 
present  joint  element.  The  problem  consists  of  a  plane  system  of  three  solia 
blocks  and  three  joint  planes.  The  system  is  restrained  at  the  bottom  and  is 
subjected  to  lateral  pressure  P  and  vertical  pressure  2P.  The  joints  are 
assumed  to  be  filled  with  nondi latent  material  having  the  properties  shown  in 
Figure  4-20.  The  blocks  are  assumed  to  have  the  isotropic  elastic  properties 
shown  in  Figure  4-20.  The  matecial  in  the  joints  is  represented  by  joint 
elements  which  are  shown  as  their  shaded  strips.  Ambiguity  with  respect  to 
nodal  displacements  would  arise  at  the  point  of  the  wedge-shaped  block  if  the 
joint  elements  were  to  extend  to  the  point.  Ambiguity  is  avoided  by  stopping 
the  elements  short  of  the  point.  This  is  a  good  approximation  to  a  physical 
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FIGURE  4-20.  GEOMETRY  OF  EXAMPLE  PROBLEM  USING  PLANE  SLIP  ELEMENT 
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system  if  the  joint  elements  exteuu  infinitesimally  close  to  the  point.  In  the 
present  example,  the  gap  is  made  laiga  to  illustrate  how  the  difficulty  is 
avoided.  The  blocks  are  represented  by  plane  strain  elements. 

Results  of  the  present  example  are  shown  in  Figures  4-21  and  4-22. 
Figure  4-21  shows  how  the  vertical  displacement  varies  as  a  function  of  hori¬ 
zontal  distance  from  the  centerline.  Each  curve  in  Figure  4-21  corresponds 

to  a  different  distance  below  the  top  surface.  Figure  4-22  shows  the  wedge  in 
its  final  pos  i  tion. 
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SECTION  5 

DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

This  section  describes  the  computer  program  in  terms  of  features 
which  are  apparent  to  the  user,  such  as  automatic  mesh  generation,  and  in 
terms  of  its  logical  structure.  The  first  part  of  the  section  descriL.es  the 
steps  which  the  user  takes  in  order  to  prepare  the  input  data.  The  mesh 
generator,  bandwidth  reducer  and  plotting  capability  are  described.  The 
motive  and  logic  diagrams  for  element  renumbering  are  also  described. 

The  second  part  of  this  section  describes  the  execution  part  of  the 
computer  program  in  terms  of  logic  diagrams.  Subroutines  are  described  which 
control  the  computation,  form  the  global  stiffness  matrix  and  compute  the  load 
vector,  form  element  stiffnesses  and  other  operations.  Also  described  is  the 
technique  of  multibuffering  whereby  data  is  transferred  between  core  and  per¬ 
ipheral  storage  units  at  the  same  time  that  computations  are  being  performed 
in  core.  This  technique  greatly  improves  the  efficiency  of  the  program  for 
problems  where  large  blocks  of  data  are  stored  on  peripheral  memory  units. 

5.1  PROCESSING  OF  INPUT  DATA 

The  sequence  in  which  input  data  are  processed  is  shown  in  Figure  5- 
There  are  three  basic  phases.  The  first  consists  of  generating  continuum 
(two-  or  three-dimensional)  and  joint  elements  and  plotting  the  results.  The 
second  phase  consists  of  adding  other  elements  (beam,  thick  shell,  bar,  etc.) 
manually,  reducing  the  bandwidth,  shuffling  the  element  numbers  so  that  they 
appear  in  the  order  of  contribution  to  the  global  stiffness  matrix,  and 
plotting  the  final  mesh.  The  third  phase  consists  of  reading  additional  data 
such  as  material  properties  and  loads. 

MESH  GENERATOR 

The  automatic  mesh  generation  scheme  incorporated  in  the  program 
requires  the  user  to  provide  a  coarse  mesh  which  the  program  refines  by 
subdivision  under  the  control  of  the  user.  The  main  goal  is  to  minimize  the 
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the  input  preparation  on  the  part  of  the  user.  The  method  uses  the  key 
diagram  concept  described  in  Reference  5“ * ;  the  remainder  of  the  present 
mesh  generation  scheme  is  new  and  original  work. 

The  subdivision  of  the  given,  coarse  mesh  to  obtain  the  final,  fine 
mesh  is  carried  out  by  subdividing  each  one  of  the  mesh  units,  as  "zones," 
w i thin  the  coarse  mesh  in  the  following  manner.  Consider  the  simple  case 
of  a  general,  two-dimensional  mesh,  in  which  the  basic  mesh  unit  or  zone  is 
a  parabolic  quadrilateral  as  shown  in  Figure  5-2.  Assume  that  the  x  and  y 
coordinates  of  the  four  corner  mesh  points,  1,  2,  3,  and  A,  as  well  as  those 
of  the  four  midpoints,  5»  6,  7i  and  b,  are  known.  The  x  and  y  coordinates 
of  an  arbitrary  point  P  within  this  zone  can  then  be  expressed  conveniently 
in  terms  of  the  coordinates  of  the  eight  points  through  the  use  of  the  local, 
curvilinear  coordinates,  £  and  n,  whose  values  range  from  -1  to  +1  on 
opposite  sides  as  indicated  in  Figure  5-2.  Thus,  one  can  write 


x 


y 


8 

D  N.x. 
i  =  1  1  1 


8 

£  Vi 

1  =  l 


(5-1) 


in  which  N.'s  are  the  so-called  "isoparametric  shape  functions"  expressible 
in  terms  of  the  curvilinear  coordinates  as  follows: 


fJ1  =  "  if  0  "  €)  0  -  n)  (€  + 

fJ2  =  "  \  0  "  S)  (1  +  n)  (5  “ 

N3  =  ^  ( 1  +  5) (1  +  n) U  +  n 

=  ^  0  +  C) (1  -  n)  (5  -  n 

N5  =  j  (1  -  00  -  nf) 

N6  =  }  (I  -  5)0  -  n2) 

j  (1  -  C2)0  +  n) 

N8  -  Y  (1  -  52)(1  +  n) 


n  +  1) 
n  +  1) 

-  1) 

-  1) 

(5-2) 
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A 

Thus,  for  a  certain  set  of  values  of  the  curvilinear  coordinates,  or  "natural 
coordinates,"  the  corresponding  Cartesian  coordinates  can  be  easily  found 
from  Equations  5~2.  Consequently,  to  subdivide  the  zone  into  a  b  x  k  mesh, 
for  example,  and  establ i sh  the  Cartesian  coordinates  of  the  twenty-one  new 
mesh  point..,  one  merely  substitutes  into  Equations  5_1  twenty-one  times, 
each  time  with  a  different  combination  of  values  of  c  and  n  which  are 
incremented  successively  by  0.5- 

This  technique  can  be  used  to  treat  each  one  of  the  zones  in  a 
given,  coarse  mesh  of  any  configuration.  The  only  restriction  that  has  to  be 
observed  is  that  the  number  of  subdivisions  between  any  two  adjoining  zones 
must  match.  This  is  to  satisfy  the  fundamental  connectivity  requirement  of 
the  finite  element  method.  To  satisfy  this  continuity  requirement  and  to 
facilitate  the  preparation  of  the  input  that  defines  the  basic  coarse  mesh, 
it  is  convenient  to  introduce  the  so-called  "key  diagram." 

A  key  diagram,  in  general,  is  a  rectangular  grid  resembling  a 
checker-board.  It  has  no  physical  dimensions.  Its  purpose  is  to  present, 
or  define,  the  connectivity  of  the  zones  in  the  coarse  mesh  and  to  facilitate 
defining  the  position  of  a  mesh  point  in  the  form  of  row  and  column  numbers 
in  the  coarse  mesh.  Another  purpose  of  the  key  diagram  is  to  help  define  the 
extent  of  subdivision  of  a  zone,  for,  to  satisfy  the  aforementioned  connec¬ 
tivity,  a  newly  introduced  mesh  line  has  to  extend  across  all  the  zones 
located  in  the  same  row  or  column  in  the  key  diagram. 

Figure  5"3  illustrates  the  generation  of  a  simple  finite  element 
mesh  representing  a  dam  and  part  of  its  foundation.  The  domain  has  been 
blocked  into  three  zones  which  are  connected  as  shown  in  the  accompanying 
key  diagram.  The  final  mesh  is  the  result  of  subdividing  the  first  and  the 
last  zone  into  two  subdivisions  and  the  second  into  three  subdivisions  in 
the  vertical  direction,  and  all  three  zones  into  three  subdivisions  in  the 
horizontal  direction. 
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The  input  to  the  program  for  this  particular  example  would  consist 
of  the  fo I  lowing: 

a.  A  statement  that  the  key  diagram,  or  the  original  coarse  mesh, 
is  3  by  1  . 

b.  The  rectangular  coordinates  of  the  given,  eight  mesh  points 
and  those  of  the  single  midpoint  that  define  the  lower  boundary 
on  a  curve. 

c.  Three  material  numbers  to  associate  the  first  two  zones  with 
concrete  and  the  third  with  the  soil  medium. 

d.  The  three  numbers,  2,  3  and  2,  that  specify  the  number  of  sub¬ 
divisions  in  the  vertical  direction,  and  the  number  2  for 

the  horizontal  direction. 

The  computer  program  would  assign  the  material  number  for  a  zone  to  all  of 
the  elements  that  are  created  in  the  zone,  and  number  the  nodal  points  and 
elements  in  the  final  mesh  column-wise,  from  the  top  to  the  bottom,  from  the 
left  toward  the  right  by  referring  to  the  key  diagram.  The  curved  mesh  lines 
in  the  third  zone  are  meant  only  to  emphasize  the  fact  that  nodal  points  are 
generated  on  a  second  degree  curve  because  of  the  nature  of  the  shape  func¬ 
tion;  the  actual  mesh  line  is  always  a  straight  line  between  two  neighboring 
mesh  points. 

In  the  previous  example,  if  the  user  is  not  satisfied  with  the 
refinement  of  the  final  mesh  or  with  that  in  any  one  of  the  zones,  he  can 
quickly  rerun  the  problem  by  modifying  the  four  integers  mentioned  in  I  tern  d 
without  having  to  change  any  of  the  data  mentioned  in  the  first  three  items 
which  define  the  given,  coarse  mesh.  This  is  perhaps  the  most  important 
feature  of  this  particular  mesh  generation  scheme. 
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Aside  from  being  used  to  define  an  edge  of  a  zone  on  a  curve,  a 
midpoint  can  also  be  used  to  imply  a  mesh  grading.  The  mesh  shown  in 
Figure  5-i<  is  the  result  of  specifying  the  three  midpoints,  Mj  ,  M^  and  , 
purposely  off-centered  and  toward  the  edge  of  the  circular  hole.  Midpoints 
M 1  and  M  are  there  merely  to  form  the  edge  of  the  hole  on  a  curve.  In 
this  particular  example,  number  of  subdivisions  of  5  and  5  in  the  vertical 
direction  (in  the  key  diagram)  and  of  16  in  the  horizontal  direction  were 
spec i f i ed . 


Another  feature  that  makes  the  scheme  versatile  and  powerful  is 
the  provision  for  the  user  to  prescribe  a  zone  as  void  merely  by  assigning 
zero  for  the  material  number  of  that  zone.  With  this  provision,  one  can 
easily  generate  a  mesh  around  a  notch  or  a  cutout.  The  mesh  around  a  tunnel 
shown  in  Figure  5“5  is  such  an  example. 

The  final  feature  to  be  mentioned  here  is  the  program's  ability  to 
join  any  two  edges  in  the  key  diagram.  This  makes  it  possible  to  produce  the 
meshes  shown  in  Figures  5~6  and  5"7-  The  user  may  prefer  the  type  of  mesh 
shown  in  Figure  5~7  over  that  shown  in  Figure  5“5  for  the  radiating  pattern 
of  the  mesh  lines  which  makes  it  easier  to  make  a  finer  mesh  around  the 
edge  of  the  tunne 1  . 

The  same  scheme  can  be  used  to  generate  a  three-dimensional  mesh. 
The  necessary  changes  are  to  replace  the  two-dimensional  shape  function  and 
key  diagram  by  their  three-dimensional  counterparts,  and  to  introduce  the 
third  Cartesian  coordinate.  The  basic  “zone"  or  "block"  in  this  case  is  a 
general,  six-sided  body  with  parabolic  edges.  It  has  eight  corner  nodes  and 
twelve  midpoints  as  shown  in  Figure  5-8>  and  the  Cartesian  coordinates  at  an 
arbitrary  point  P  within  the  block  having  the  natural  coordinates  £ ,  n 
and  <;  are  expressible  in  terms  of  the  Cartesian  coordinates  of  the  twenty 
controlling  points  as  follows: 
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x  = 

y  = 


where,  for  corner  points, 

N.  =  i  (1  +  c.c)(l  +  n.n)  (1  +  c . c) (c ; C  +  n.n  +  -  2)  (5-4 

10  I  I  I  I  t  I 

and  for  midpoints, 

N.  =  ^  (1  -  ?2)(1  +  n  j  n)  ( 1  +  C j  c)  for  points  with  ?.  =  0 

i  2 

Nj  =  -Jj-  (1  +  5.5)(1  -  n  )(1  +  SjO  for  points  with  n j  =  0 

i  2 

Nj  =  ^  (1  +  5.00  +  n-n)  0  "  5  )  for  po.ints  with  c.  =  0 

As  before,  zero  can  be  used  as  material  number  to  imply  void  or  cavity,  and 
a  midpoint  can  be  so  specified  as  to  imply  a  curved  edge  as  well  as  a  mesh 
grading  along  the  edge.  Instead  of  joining  two  edges,  one  can  now  join  any 
two  interfaces  as  long  as  the  meshes  on  them  are  compatible.  Some  three- 
dimensional  meshes  automatically  generated  in  this  manner  are  illustrated 
in  Figures  5-9  and  5~10. 
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BANDWIDTH  REDUCER 

The  computer  time  required  to  perform  a  computation  is  approximately 
proportional  to  the  square  of  the  bandwidth.  Thus  it  is  important  to  aim 
for  a  minimum  bandwidth  when  numbering  nodal  points.  However,  when  the  finite 
element  mesh  is  generated  automatically  by  the  technique  described  above, 
no  attention  whatever  is  paid  to  minimizing  bandwidth.  The  user  can  influence 
the  bandwidth  to  some  extent  by  judicious  choice  of  key  diagram.  Neverthe¬ 
less,  it  is  very  likely  that  the  bandwidth  so  generated  will  not  be  optimum. 
The  situation  will  often  be  made  worse  when  elements  are  added  manually. 

To  help  alleviate  this  difficulty  and  thus  encourage  the  user  to  use  the 
automatic  mesh  generator,  a  bandwidth  reducer  is  included  in  the  present 
computer  program. 

The  function  of  the  bandwidth  reducer  is  shown  in  Figure  5 “ 1 1 .  The 
INPUT  NODAL  POINT  configuration  is  typical  of  the  mesh  which  would  be  gen¬ 
erated  automatically  for  which  the  maximum  difference  in  nodal  point  numbers 
for  any  element  is  31.  This  configuration  was  submitted  to  the  bandwidth 
reducer.  The  result,  after  2  sec  of  computation  time  by  the  bandwidth 
reducer  is  the  REDUCED  NODAL  POINT  configuration,  for  which  the  maximum 
difference  in  node  numbers  is  6.  That  this  configuration  is  not  optimum  is 
shown  by  the  IDEAL  configuration  for  which  the  maximum  difference  is  5-  Thus, 
the  technique  is  not  optimum  because  it  does  not  always  converge  to  the 
minimum  bandwidth  in  the  computer  time  which  the  user  has  selected.  Also, 
the  technique  operates  on  nodal  point  numbers  rather  than  on  degree  of 
freedom  numbers.  Thus,  when  several  types  of  elements  are  mixed,  the  con¬ 
figuration  corresponding  to  minimum  difference  in  node  numbers  does  not 
necessarily  correspond  to  minimum  bandwidth. 

Logic  diagrams  for  the  bandwidth  reducer  are  shown  in  Figures  5-12 

and  5-13. 

ELEMENT  NUMBERING 

In  addition  to  renumbering  nodal  points  to  reduce  the  bandwidth  of 
the  global  stiffness  matrix,  the  elements  are  numbered  in  the  order  of  their 
contribution  to  the  global  stiffness  matrix.  In  this  way,  the  strain/ 


R-7215-1-2701 


INPUT  NODAL  POINT  CONFIGURATION 
MAXIMUM  NODAL  POINT  DIFFERENCE  =  31 


REDUCED  NODAL  POINT  CONFIGURATION 
MAXIMUM  NODAL  POINT  DIFFERENCE  -  6 
SOLUTION  TIME  =  2  SECONDS 


IDEAL  CONFIGURATION 

MAXIMUM  NODAL  POINT  DIFFERENCE  =  5 


AA  1(681 


FIGURE  5-11.  EXAMPLE  OF  BANDWIDTH  REDUCER 
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FIGURE  5-12.  PROCEDURE  USED  IN  BANDWIDTH  REDUCER  (REFERENCE  5-2) 
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stiffness  matrixes  can  be  retrieved  from  peripheral  storage  in  the  least 
time.  iince  the  el  ment  data  are  used  in  the  basic  operations  of  forming 
the  global  stiffness  matrix  and  effective  load  vector,  the  efficiency  gained 
by  performing  these  complicated  operations  is  essential. 

Tlie  element  numbers  assigned  by  the  mesh  generator  and  the  user 
are  revised  such  that 

f  k 

mi n  (LM ( I ) )  for  the  n  element  < 
mi  n  ( LM ( I ) )  for  the  n  +  1^  element 


where 

L-M(l)  =  Array  of  degree  of  freedom  members  for  the  element 

Element  data  are  stored  sequentially  on  peripheral  storage  units  such  that 
the  data  for  Element  1  is  at  the  head  of  the  unit,  and  then  in  sequence 

Element  Number  =  1  <  . ..  n,  n  +  1  ...  <  NUMEL 


where 

NUMEL  =  Total  number  of  elements 

In  order  to  modify  the  effective  load  vector  F,  to  obtain  the 
matrix  C  of  s t ress/s tra i n  coefficients  and  to  obtain  element  stresses  and 
strains,  an  array  containing  the  previous  loads  F,  the  displacements  u  and 
the  incremental  displacements  du  is  brought  into  core  as  shown  in  Figure  5-15. 
Also  required  are  data  stored  on  the  element  data  tape  such  as  the  strain/ 
displacement  transformation  matrixes  and  LM  arrays  for  the  elements  of 
interest.  As  shown  in  Figure  F  i*),  all  degrees  of  freedom  to  which  Element  n 
contributes  are  found  in  a  sequence  of  F,  j,  du  array  one  bandwidth  (MBAND) 
long.  Thus  it  is  possible  to  process  to  F,  u,  du  array  by  reading  the 
element  data  tape  only  once. 
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In  order  to  assemble  the  global  stiffness  matrix,  which  is  stored 
in  blocks,  Block  nn  is  brought  into  core  as  shown  in  Figure  5-15.  At  the 
same  time,  data  from  Element  n  is  made  available  from  the  element  data  tape. 
Element  n  fits  entirely  into  Block  nn .  Elements  n,  nn  +  1,  etc.  are 
processed  so  long  as 

DFE1  >  DFI<2 

and 

dfe2  <  dfk2 

In  this  event,  element  stiffness  is  directly  added  to  the  appropriate  rows  and 
columns  of  Block  nn.  If 

dfe]  >  dfi<2 

and 

dfe2  >  dfk2 

it  signifies  that  Element  n  contributes  to  Block  nn  and  to  Block  nn  +  1  . 

In  this  event,  those  contributions  which  can  be  made  to  Block  nn  are  made, 
and  the  stiffnesses  and  LM  arrays  of  those  elements  which  contribute  also 
to  Block  nn  +  1  are  written  temporarily  on  storage  unit  T1  for  subsequent 
insertion  into  Block  nn  +  1 .  Elements  are  processed  and  their  stiffnesses 
are  added  to  Block  nn  and  or  accumulated  on  T1  until 

DFE]  >  DFK2 

This  signifies  that  Block  nn  of  the  global  stiffness  matrix  is  prepared  except 
for  data  stored  on  unit  TO,  which  contains  element  data  overlapping 
Blocks  nn  -  1  and  nn.  This  data  is  read  into  core  and  that  part  for  which 

DFE 1  <  DFK] 
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is  added  to  Block  nn.  Now  Block  nn  is  completely  prepared  and  is  moved  out 
of  core  while  Block  nn  +  1  is  brought  into  com.  Block  nn  +  1  is  filled  in 
the  same  way  as  Blocks  nn  -  I  and  nn,  which  is  by  passing  sequentially  through 
the  element  data  array.  When  all  direct  conn  i  but  ion.,  to  Block  nn  +  1  have 
been  made  and  all  data  which  contributes  to  Block  nn  +  2  has  been  stored  on 
unit  T2,  the  data  stored  on  unit  T1  is  read  into  core  and  added  to 
Block  nn  +  1 .  This  block  is  now  completely  prepared  and  is  transferred  to 
peripheral  storage. 

SIMULATION  OF  THE  CONSTRUCTION  AND  EXCAVATION  SEQUENCE 

The  construction  and  excavation  sequence  will  be  represented  by 
modification  of  the  material  properties  of  the  elements  involved.  Initially 
elements  will  be  assigned  to  the  regions  to  be  constructed  or  excavated  later. 

A  flag  will  indicate  the  time  of  the  addition  or  subtraction  of  each  element. 

In  case  of  excavation,  the  elements  of  that  region  will  have  the  appropriate 
material  properties  and  will  contribute  to  the  global  stiffness  of  the  system 
up  to  the  indicated  time,  beyond  which  the  contribution  of  these  elements  to 
the  global  stiffness  of  the  system  will  be  zero.  The  reverse  of  this  procedure 
will  be  applied  to  the  elements  of  the  regions  to  be  constructed  later.  All 
these  operations  will  be  performed  in  the  element  package. 

SELF-LOADING 

The  initial  state  of  the  system  will  be  assumed  to  be  stress  free. 
Then  the  dead  load,  if  any,  will  be  applied  in  a  specified  number  of  increments 
prior  to  the  application  of  the  external  load.  It  is  necessary  to  apply  the 
self-loading  incrementally  due  to  the  fact  that  the  system  is,  in  general, 
non  1 i nea  r . 
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5.2  EXECUTION  PHASE  OF  PROGRAM 

The  execution  phase  of  the  program  is  summarized  in  Figures  5-16 
through  5-21.  The  controlling  program  is  BMCALC,  Figure  5- 1 6 .  The  first 
main  operation  is  to  form  the  effec.ive  load  vector  and  global  stiffness 
matrix,  which  is  performed  by  Subroutine  KFORM ,  Figure  5-17.  The  assembly 
of  the  global  stiffness  from  t lie  element  stiffnesses  is  performed  by  Sub¬ 
routine  BSTIF,  Figure  5-1 8 .  Subroutines  FWRT ,  TDRUM  and  FILLFU  transfer  data 
from  peripheral  storage  to  core. 

5.2.1  MULTIBUFFERING  TECHNIQUE 

Mu  1 t i bu f fer i ng  is  a  technique  whereby  central  processor  wait  time 
for  all  binary  Reac /Write  operations  involving  the  peripheral  storage  of  data 
is  minimized.  (Formatted  I/O  operations  such  as  card  reading,  punching,  and 
printing  are  not  included  in  this  discussion.)  The  reason  for  using  multibuf¬ 
fering  techniques  is  that  data  moves  faster  between  core  locations  than  between 
peripheral  storage  and  core. 

The  problem  which  mu  1 1 i buf f er i ng  overcomes  is  the  standard  I/O 
feature  of  higher  level  programming  languages,  such  as  FORTRAN,  which  requires 
that  when  an  I/O  operation  (READ  or  WRITE)  is  initiated,  computation  ceases 
until  the  I/O  operation  is  completed.  This  feature  assures  the  user  that  data 
he  may  wish  to  use  is  in  core  before  he  tries  to  use  it.  The  amount  of  time 
the  program  must  wait  for  completion  of  1/0  operations  depends  on  (1)  the 
access  time  and  (2)  the  data  transfer  time.  These  times  depend  on  the  type 
of  peripheral  device  being  used  and  the  amount  of  data  to  be  transferred. 

Thus,  as  the  amount  of  peripheral  storage  increases,  the  time  spent  waiting 
for  completion  of  I/O  operations  increases  and  for  large  volumes  of  data 
the  I/O  time  may  control  overall  run  times.  Mu  1 1 i buf fer i nq  minimizes  the 
wait  time  of  these  I/O  operations  by  allowing  computations  to  proceed  at  the 
same  time  data  transfer  from  peripheral  storage  is  occurring.  This  requires 
standard  FORTRAN  I/O  operations  to  be  replaced.  This  is  possible  on  most 
large  scale  scientific  computers  by  the  use  of  special  machine-dependent 
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FIGURE  5-16.  PROGRAM  BMCALC- -CONTROLS  MAIN  OPERATIONS  OF  THE 
COMPUTATION  SECTION 
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FIGURE  3-17.  SUBROUTINE  KFORM--CALLS  SUBROUTINES  TO  COMPUTE  THE  LOAD 
VECTOR  AND  GLOBAL  STIFFNESS 
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FIGURE  5-17.  (CONTINUED) 
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BEAM  ELEMENTS  THICK  SHELL  ELEMENTS  l-D,  2-D,  AND  3-D  ELEMENTS 

[k]  IS  WITH  ELEMENT  DATA  [k]  IS  ON  ELEMENT  DATA 

FILE  (S  KELOKDS) 


FIGURE  5-18.  SUBROUTINE  BST I F--COMPUTER  ELEMENT  STIFFNESS  [k]  AND  ADDS 
TO  GLOBAL  [k] 
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FIGURE  5-19.  SUBROUTINE  FWRT-MOVES  DATA  FROM  LIVE  LOAD  VECTOR  (FWORK) 
TO  F,  u,  du  OUTPUT  BUFFER  AND  TO  u  OUTPUT  FILE 


R- 7215-1-2701 


AA4666 


FIGURE  5-20. 


SUBROUTINE  TDRUM--ADDS  IN  ELEMENT  [k]  WHICH  OVERFLOWED  FROM 
PREVIOUS  BLOCK  OF  [k] 
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statements  or  subroutines.  These  interface  routines  are  generally  different 
on  each  machine  and,  unless  care  is  exercised,  a  program  may  become  machine 
dependent.  This  danger  is  avoided  in  the  present  application  of  multibuffer¬ 
ing  by  isolating  all  I/O  statements  in  one  subroutine  which  may  easily  be 
modified  when  moving  the  program  to  a  different  machine. 

A  general  mu  1 t i buf f er i nq  scheme  to  perform  computat'ons  on  a  set  of 
data  stored  on  peripheral  devic'’*  is  shown  in  Figure  5"22.  This  scheme, 
which  is  incorporated  in  the  present  computer  procram,  requires  either  the 
amount  of  main  storage  used  for  buffers  to  he  increased  or  the  s i ze  of  each 
buffer  to  be  decreased  relative  to  the  buffer  size  that  might  be  used  with 
standard  FORTRAN  I/O  procedures.  Since  most  of  the  main  storage  is  alreadv 
used  in  determining  buffer  sizes,  the  latter  alternative  is  employed.  fxtensive 
testing  and  previous  experience  have  shown  that,  although  the  number  of  I/O 
operations  increases,  mu  I t i buf f er  i  ng  results  in  substantial  overall  reduction 
in  computer  run  time  for  a  given  problem.  Savings  increase  with  the  problem 
size  and  thus  the  volume  of  data  on  peripheral  storage  increases. 

5.2.2  BAND  SOLVER 

The  solution  of  larqe  structural  systems  requires  he  solution  of 
a  set  of  linear  simultaneous  equations  of  the  form 

i F}  =  [K]  fu)  (5-5) 


where 

{ F }  is  a  vector  of  applied  loads 

{ u }  is  a  vector  of  unknown  displacements  or,  as  in  the  present 
case,  displacement  increments 

f K }  is  the  global  stiffness  matrix.  In  the  present  case  it  is 
banded  and  positive  definite 
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FIGURE  5-22.  TYPICAL  OPERATION  OF  MULT  I  BUFFER  I  NG  TECHNIOLIF 
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Many  methods  of  obtaining  a  so  1  u  ion  to  fq/ation  5"5  are  available.  A 
frequently  used  method  is  the  Cholt.ki  Decomposition  Method.  Defining 

IK]  =  [LI  [D]  [L]T  (5-6) 


where 


[L] 

i  s 

[D] 

i  s 

i  tut 

i  nq 

IF/ 

_ 

i  x 


and  defining 

=  [D]  [ L] T  f u } 

Equation  5-7  becomes 

(  Ft  =  [L]  { Z } 

There  are  many  algorithms  for  solving  Equations  5-8  and  5-9, 
alg<  1  ithms  used  in  this  code  are  (Reference  5~^0 
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where 

n  =  Total  number  of  equations 

i  ,j  =  Row  and  column  indices 

The  use  of  Equations  5“ 1 0  and  5-11  to  obtain  the  [L]  and  [ D) 
matrices  in  the  most  general  case  of  a  block-by-block  solution  is  illustrated 
in  Figures  5 - 2 3  through  5~ 2 5 •  Figure  5- 2 3  shows  a  typical  banded  stiffness 
matrix  with  the  arrangement  of  the  terms  in  each  core  block.  Although  shown 
as  a  two-dimensional  array,  the  actual  locations  of  storage  may  be  consecutive. 
Figure  5-2*4  shows  a  method  of  assigning  core  storage  to  allow  the  double 
buffering  scheme  of  Section  5-2.1  to., be  used.  The  indicated  Scratch  Core 
area  may  be  of  any  size  (it  must  be  at  least  four  words  in  length)  and  is 
used  as  a  buffer  area  to  save  columns  of  the  reduced  stiffness  matrix  needed  for 
for  reduction,  i.e.,  decomposition  of  future  stiffness  blocks.  All  I/O  opera¬ 
tions  between  the  scratch  area  of  core  and  peripheral  devices  use  the  technique 
of  Section  5-2.1.  Thus  storage  and  retrieval  of  intermediate  data  required  for 
the  computation  of  [L]  and  [D]  may  in  the  optimum  case  be  performed  at 
core  speed.  Only  a  single  output  buffer  is  shown  in  Figure  5-2*4.  This  presents 
no  contradiction  to  the  double  buffering  scheme.  Figure  5~ 2 5  shows  that 
between  initiation  of  the  write  and  storage  of  new  data  in  the  output  buffer 
many  calculations  are  performed.  Tests  have  shown  that  this  calculation  time 
is  greater  than  the  data  transfer  time  thus  allowing  a  single  output  buffer 
to  be  used.  Figure  5  —  2 5  shows  the  sequence  of  operations  in  decomposing  the 
stiffness  matrix.  The  generation  of  { Z }  (from  Equation  5 ~ 1 2 )  may  be  per¬ 
formed  at  the  same  time  [L]  and  [D]  are  generated. 

After  completing  decomposition  of  the  stiffness  matrix  and  generating 
the  { Z }  matrix,  the  displacements  {u}  are  computed  by  Equation  5-13. 

Figure  5-26  shows  the  sequence  of  operations  in  computing  the  { u }  matrix. 
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FIGURE  5- 


INPUT  STIFFNESS  BUFFER  tf  1 


INPUT  STIFFNESS  BUFFER  HI 


OUTPUT  STIFFNESS  (REDUCED) 
BUFFER 


SCRATCH  CORE  AREA 


CORE  BUFFERS  FOR  STIFFNESS  MATRIX  DECOMPOSITION 
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FIGURE  5-25- 


SCHEMATIC  DIAGRAM  OF  STIFFNESS  MATRIX  DECOMPOSITION 
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FIGURE  5-26.  SCHEMATIC  DIAGRAM  OF  SOLUTION  VECTOR  EVALUATION 
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SECTION  6 

COMPARISON  WITH  CLOSED-FORM  ANALYT  C  SOLUTIONS 


6. 1  SAMPLE  PROBLEMS 

To  investigate  the  numerical  accuracy  of  the  present  computer  program 
and  to  determine  the  computer  time  required  to  solve  problems  of  various  size 
and  complexity,  several  example  problems  have  been  formulated  and  their 
numerical  solutions  compared  with  closed-form,  analytic  solutions.  One-,  two-, 
and  three-dimensional  linear  elastic  solutions  are  considered  as  well  as  one- 
and  two-dimensional  elastic/plastic  and  two-dimensional  visco-elastic  solu¬ 
tions.  These  are  listed  in  Table  6-1. 

Problem  1--Stress  Around  a  Circular  Ho  1 e 

The  geometry  of  Problem  No.  1  is  illustrated  in  Figure  6-1 .  The 

finite  element  mesh  is  illustrated  in  Figure  6-2.  The  solution  is  shown  in 

Figure  6-3  in  te--ms  of  principal  stresses  at  0  =  5.7  and  ^2.1°.  It  does  not 
depend  on  the  material  properties  of  the  plate. 

Problem  2--Stress  and  Displacement  in  an  Elastic 
Thick-Wal led  Cylinder  Under  Internal  Pressure 

The  geometry  of  Problem  No.  2  is  illustrated  in  Figure  6-^4.  The 

finite  element  mesh  is  also  shown.  The  solution  is  shown  in  Figure  6-5  in 

terms  of  radial  and  tangential  stresses. 


Problem  3~-Stress  in  an  Elastic,  Perfectly  Plastic 
Thick-Walled  Cylinder  Under  Internal  Pressure 

The  geometry  of  Problem  No.  3  is  the  same  as  that  of  Problem  No.  2 
and  is  shown  in  Figure  6-A.  For  this  problem,  an  additional  material  property, 
the  Mises  yield  criterion,  is  specified  as  follows: 


f  =  -  a,  <  0 


2  “1  - 


(6-1) 


The  analytic  solution  is  shown  in  Figure  6-6. 
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TABLE  6-1.  PROBLEMS  SOLVED  BY  FINITE  ELEMENT 
AND  CLOSED-FORM  METHODS 


Stress  around  a 
ci rcul ar  hole 

Two-di mens i onal 
(p 1 ane) 

E 1  as  t  i  c 

Stress  in  a  thick 
walled  cy 1 i nder 
under  internal 
pressure 

One-d i mens iona 1 
(axi symmetri c) 

Elastic 

Stress  in  a  thick 
walled  cy 1 i nder 
under  internal 
pressure 

One-d i mens i ona 1 
(axi symmetri c) 

Elastic,  perfectly 
pi  as  ti c 

Stress  in  a  rein¬ 
forced  thick 
wal led  cyl i nder 

One-d i mens iona 1 
(axi symmetri c) 

V i scoel as t i c 

Stress  concentra¬ 
tion  around  a 
cylindrical  hole 
i n  a  semi  infinite 
body 

Three- 

dimensional 

Elastic 

Closed-Form 
Sol ut  ion 

Reference  6-1 

Reference  6-2 

Reference  6-2 

Reference  6-3 

Reference  6-^4 
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SQUARE  PLATE  (  20  IN.  X  20  IN.  ) 


DISTANCE  (  IN) 


FIGURE  6-2.  FINITE  ELEMENT  MESH  FOR  STRESS  CONCENTRATION  AROUND 
CIRCULAR  HOLE 


NORMALIZED  PRINCIPAL  STRESS 
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a)  ACTUAL  GEOMETRY  (b)  FINITE  ELEMENT  REPRESENTATION 
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Problem  4--Stress  in  a  Visco-Elastic,  Reinforced 
Cylinder  Under  Internal  Pressure 

The  geometry  of  Problem  No.  4,  shown  in  Figure  6—  7 >  is  similar  to 
that  of  Problems  2  and  3.  The  main  difference  is  a  steel  reinforcing  ring 
around  the  outer  circumference.  The  material  of  the  cylinder  is  assumed  to  be 
governed  by  a  Maxwell -type  law  as  follows: 


j.j  =  2G(e.j)  exp  ( - 1/ B ) 


(6-2) 


where 

c!  .  =  Component  of  deviatoric  strain  tensor 

G  =  Elastic  shear  modulus 

t  =  Time  in  units  of  B  =  n/G  where  n  -  viscoelastic  parameter) 

Volumetric  deformation  of  the  cylinder  and  all  deformations  of  the  reinforcing 
ring  are  assumed  to  be  linearly  elastic  inviscid.  The  variation  of  radial  and 
circumferential  stresses  are  shown  as  functions  of  radius  and  time  in  Figure  6-8 

Problem  5~~Three-D imens !ona 1  Stress  Concentration  Around  a 
Cylindrical  Hole  it:  a  Semiinfinite  Elastic  Body 

The  geometry  of  Problem  No.  5  is  illustrated  in  Figure  6-9.  The 
stress  distribution  around  the  hole  near  the  stress  free  face  (X^  -  X2  plane) 
is  appropriate  to  plane  stress,  while  in  the  interior  there  is  axial  stress 
along  the  axis  of  the  hole.  The  finite  element  mesh  is  shown  in  Figure  6-10. 

The  loading  condition  selected  for  this  example  is  uniaxial  stress  parallel 
to  the  X-axis.  Thus  the  faces  parallel  to  the  X2~p  plane  are  stress  free  as 
is  one  face  parallel  to  the  X2-X^  plane.  The  finite  element  solution  is 
compared  with  the  analytic  solution  in  Figure  6-11. 


VISCO  ELASTIC  MATERIAL 


PROBLEM  '4— INFINITELY  LONG  REINFORCED  THICK  VISCOELASTIC  CYLINDER 
SUBJECTED  TO  INTERNAL  PRESSURE 


FIGURE  6-o.  PROBLEM  4--RE I NFORCED  VISCOELASTIC  CYLINDER  SUBJECTED  TO  INTERNAL 
PRESSURE  (REFERENCE  6-4) 


PROBLEM  5--FINITE  ELEMENT  MESH  (ONLY  FIRST  THREE  U»V*RS  ARE 
SHOWN  FOR  CLARITY.  COMPLETE  MESH  CONTAINS  EIGHT 
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A 

6.2  COMPUTING  TIME  REQUIRED  FOR  SOLUTION 

The  computing  time  required  to  solve  Problem  5  and  several  small 
problems  on  a  Univac  1108  with  65,000  words  of  core  is  shown  in  Table  6-2. 
More  data  on  time  to  solve  various  sizes  of  problems  will  be  gathered  during 
the  second  phase  of  the  contract. 
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APPENDIX  A 

LOGIC  DIAGRAMS 
FOR  MATERIAL  PROPERTIES 

Logic  diagrams  for  subroutines  in  the  material  property  package 
are  shown  below.  Subroutine  CONECT  connects  the  package  to  the  main  program. 
Subroutine  ELPL  controls  all  the  other  subroutines. 


A- 1 


FIGURE  A-1.  C ON ECT- -SUBROUTINE  CONECT  LINKS  MAIN  LINE  PROGRAM  TO 
MATERIAL  PACKAGE  ELPL 

A-2 
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FIGURE  A-2 .  SUBROUTINE  ELPL  CONTROLS  MATERIAL  PROPERTY  SUBROUTINES,  TRANSFERS 
AXES  FOR  ANISOTROPIC  MATERIALS,  PERFORMS  TESTS  FOR  STRAIN 
SPLITTING  AND  ITERATION. 


A-3 


R-721 5- 1-2701 


SPLITTING  CHECK 


mm  •  0 

Qt  -  rldoi 


Q2  -  Ql 

CONVQ  -  C0NV*Q2 
NUH  ■  HUH  ♦  1 


CALL  TSTYLO  FOR  FACTO* 


1 1- FAC TOAj  5  CQNV 


NUH  *  NO,  OF  ITFRATIONS 


Lno  (splitting 

S  REQUIRED) 


CALL  INCSTR  FOR  ELASTIC  ANO 
PLASTIC  STRAIN  INCREHENTS 


CALL  ELFUN  FOR  ELASTIC  CONSTANTS 


CALL  ELAST  FOR  Cc. 
C!  -  CEL 


CALL  YLOFUN  FOR  Cp 


OLO 

C2  -  (C  ♦  C1)/2 


CALL  ANISOT  FOR  PLANE  STRESS 
CORRECTIONS 


CALL  COHPUT  FOR  NEW  ESTIMATE 
OF  STRESSES 


QQ  -  |fe 


r  NO  (SPLITTING  REQUIREO) 


UPOATE  ELASTIC  STRAIN  VARIABLES, 
ETC.  UMAX,  VMAX,  do 


CALL  INCSTR  FOR  ELASTIC  ANO 
PLASTIC  STRAIN  INCREMENTS 


Ql  •  Z |do| 


IQ*  *  02 1  •  CONVQ 


ITERATION  LOOP 


CALL  ANISOT  TO  TRANSFORM  STRESS  I  NO 
TO  GLOBAL  COOROINATE  I 


CALL  TSTYLO  FOR  FACTOR  WHICH  IS 
TO  BE  PRINTEO  FOR  THIS  ELEMENT 


FIGURE  A-2  (CONTINUED) 


ISOTROPIC  MATERIAL 


FIGURE  A-3.  SUBROUTINE  ELAST--ELAST  FORMULATES  C-MATRIX  USING  COEFFICIENTS  FOR 
EITHER  ISOTROPIC  OR  ANISOTROPIC  MATERIALS  GENERATED  BY  ELFUN 
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FIGURE  A- 


START 


SUBROUTINE  PLAST— ' THIS  SUBROUTINE  MODIFIES  C  MATRIX 
GENERATED  BY  SUBROUTINE  ELAST  IN  ORDER  10  ACCOUNT  FOR 
PLASTICITY.  THE  BASIC  QUANTITIES  FOR  THIS  MODIFICATION 
ARE  DERIVATIVES  OF  THE  YIELD  FUNCTION  WITH  RESPECT  TO 
THEIR  ARGUMENTS  AND  ARE  COMPUTED  IN  YLDFUN. 


I\-C 


START 


St  T  I  RIG  -  0 

IRtG  -  0  IMPLItS  MATERIAL 
IS  INITIALLY  ASSUMEO  ELASTIC. 
IF  YIELOING  OCCURS  ON  PLASTIC 
SURFACE,  IREG  IS  SET  TO  I. 

If  YIELOING  OCCURS  ON  CAP, 

IREG  IS  SET  TO  2. _ 


COMPUTE  STRESS  INVARIANTS 
CJ1  ANO  CJ7 

(CJ1--IST  STRESS  INVARIANT; 
{CJ2)2--?NO  INVARIANT  OF  STRLS' 
OEVIATOR,  YCOF(I)  TO  YCOF (9) - 
COEFFICIENTS  OE FINING  ORIENTA¬ 
TION  ANO  OtGREt  OF  ANISOTROPY. 
THESi  ARE  ALL  EQUAL  TO  \  FOR 
I SOTROPIC  MODEL.) _ 


HYPI  •  ft  Of  (70) 

ITTPI--A  HAG  TO  eiiVIHCtMSH 
Dirff«4«T  TtPtS  0#  PlKrUTlY 
PUWtC  Villi  '.U*/AIIS 


SUBROUTINE  TSTYLO-- COMP ARES  CURRENT  STATE  OF  STRESS  WITH  YIELD 
CRITERIA  TO  DETERMINE  WHETHER  YIELDING  OCCURS 
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START 


RETURN 


FIGURE  A-6.  SUBROUTINE  YLDFUN--TH I S  SUBROUTINE  CALCULATES  THE  DERIVATIVES 
OF  THE  YIELD  FUNCTION  WITH  RESPECT  TO  THE  STRESS  COMPONENTS. 

THE  YIELD  FUNCTION  WHOSE  DERIVATIVES  ARE  COMPUTED  IS  DETERMINED 
BY  THE  INDICES  IREG  (SET  IN  (TSTYLD)  AND  ITYPE  (SET  BY  INPUT)) 
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IS  MATERIAL  I  ANISOTROPIC  MATERIAL 

iSOTRCPlC?  ^00"^  H  VARIABLE  E, 
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CALCULATE  STRESS 
RELAXATION  OUE  TO 
VISCOPLASTICITY 


RETURN 


I  CHECK  =  2 
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FIGURE  A-9. 


SUBROUTINE  VPLAST--TH I S  SUBROUTINE  IS  USED  FOR  A  VISCOPLASTIC 
MATERIAL.  IT  CALCULATES  THE  STATIC  YIELD  CRITERION,  DETERMINES 
WHETHER  YIELDING  OCCURS  AND,  IF  YIELDING  OCCURS,  MODIFIES  THE 
STRESSES  TO  ACCOUNT  FOR  VISCOPLASTIC  STRESS  RELAXATION. 
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START 


FIGURE  A- 1 0 .  SUBROUTINE  COMPUT--COMPUT  CALCULATES  STRESS  INCREMENTS  AND  ADDS 

THEM  TO  OLD  STRESSES  TO  OBTAIN  NEW  STRESSES.  STRAINS  CORRESPONDING 
TO  ZERO  STRESSES  (PLANE  STRESS,  UNIAXIAL  STRESS) ARE  COMPUTED. 
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